diff --git a/.gitignore b/.gitignore index 194f112..b416265 100644 --- a/.gitignore +++ b/.gitignore @@ -79,3 +79,6 @@ bower_components/ .grunt/ src/vendor/ dist/ + +# CMake artifact generation directory +/cmake_build diff --git a/CMakeLists.txt b/CMakeLists.txt index db36567..634ee80 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -55,6 +55,13 @@ set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -I/usr/X11/include") # add the executable #message("CXXLIBS" ${CMAKE_CXX_FLAGS}) +# ----------------------------------------------------------------------------- +# Find the Boost library and configure usage +set(Boost_USE_STATIC_LIBS OFF) +set(Boost_USE_MULTITHREADED ON) +set(Boost_USE_STATIC_RUNTIME OFF) +find_package(Boost 1.79.0 REQUIRED COMPONENTS serialization) + IF((NOT ${NGEN}) AND (NOT ${STANDALONE})) message("Pseudo-framework build") @@ -62,14 +69,16 @@ ENDIF((NOT ${NGEN}) AND (NOT ${STANDALONE})) if(STANDALONE) add_executable(${exe_name} ./src/bmi_main_lgar.cxx ./src/bmi_lgar.cxx ./src/lgar.cxx ./src/soil_funcs.cxx - ./src/linked_list.cxx ./src/mem_funcs.cxx ./src/util_funcs.cxx ./src/aet.cxx + ./src/linked_list.cxx ./src/mem_funcs.cxx ./src/util_funcs.cxx ./src/aet.cxx ./src/Logger.cpp ./include/Logger.hpp ./giuh/giuh.h ./giuh/giuh.c) target_link_libraries(${exe_name} PRIVATE m) + target_link_libraries(${exe_name} PRIVATE Boost::serialization) elseif(UNITTEST) add_executable(${exe_name} ./tests/main_unit_test_bmi.cxx ./src/bmi_lgar.cxx ./src/lgar.cxx ./src/soil_funcs.cxx - ./src/linked_list.cxx ./src/mem_funcs.cxx ./src/util_funcs.cxx ./src/aet.cxx ./giuh/giuh.h + ./src/linked_list.cxx ./src/mem_funcs.cxx ./src/util_funcs.cxx ./src/aet.cxx ./src/Logger.cpp ./include/Logger.hpp ./giuh/giuh.h ./giuh/giuh.c) target_link_libraries(${exe_name} PRIVATE m) + target_link_libraries(${exe_name} PRIVATE Boost::serialization) endif() @@ -82,18 +91,21 @@ add_compile_definitions(BMI_ACTIVE) if(WIN32) add_library(lasambmi SHARED src/bmi_lgar.cxx src/lgar.cxx ./src/soil_funcs.cxx ./src/linked_list.cxx ./src/mem_funcs.cxx - ./src/util_funcs.cxx ./src/aet.cxx ./giuh/giuh.c include/all.hxx ./giuh/giuh.h) + ./src/util_funcs.cxx ./src/aet.cxx ./src/Logger.cpp ./include/Logger.hpp ./giuh/giuh.c include/all.hxx ./giuh/giuh.h) else() add_library(lasambmi SHARED src/bmi_lgar.cxx src/lgar.cxx ./src/soil_funcs.cxx ./src/linked_list.cxx ./src/mem_funcs.cxx - ./src/util_funcs.cxx ./src/aet.cxx ./giuh/giuh.c include/all.hxx ./giuh/giuh.h) + ./src/util_funcs.cxx ./src/aet.cxx ./src/Logger.cpp ./include/Logger.hpp ./giuh/giuh.c include/all.hxx ./giuh/giuh.h) endif() + target_include_directories(lasambmi PRIVATE include) set_target_properties(lasambmi PROPERTIES VERSION ${PROJECT_VERSION}) set_target_properties(lasambmi PROPERTIES PUBLIC_HEADER ./include/bmi_lgar.hxx) +target_link_libraries(lasambmi PRIVATE Boost::serialization) + include(GNUInstallDirs) install(TARGETS lasambmi diff --git a/LICENSE b/LICENSE index 6c62c71..ae4bebc 100644 --- a/LICENSE +++ b/LICENSE @@ -1,3 +1,26 @@ +Copyright 2025 Raytheon Company + +Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation +and/or other materials provided with the distribution. + +Licensed under: https://opensource.org/license/bsd-2-clause + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +All rights reserved. Based on Government sponsored work under contract GS-35F-204GA. + +----------------- + +“Software code created by U.S. Government employees is not subject to copyright +in the United States (17 U.S.C. §105). The United States/Department of Commerce +reserve all rights to seek and obtain copyright protection in countries other +than the United States for Software authored in its entirety by the Department +of Commerce. To this end, the Department of Commerce hereby grants to Recipient +a royalty-free, nonexclusive license to use, copy, and create derivative works +of the Software outside of the United States.” + Apache License Version 2.0, January 2004 http://www.apache.org/licenses/ diff --git a/include/Logger.hpp b/include/Logger.hpp new file mode 100644 index 0000000..c884c3e --- /dev/null +++ b/include/Logger.hpp @@ -0,0 +1,66 @@ +#ifndef LOGGER_HPP +#define LOGGER_HPP + +#include +#include +#include +#include +#include +#include +#include + +#define LOG Logger::Log + +enum class LogLevel { + NONE = 0, + DEBUG = 1, + INFO = 2, + WARNING = 3, + SEVERE = 4, + FATAL = 5, +}; + +/** + * Logger Class Used to Output Details of Current Application Flow + All methods and variables are static so instantiating an object is unnecessary. + */ +class Logger { + public: + // Methods + static void Log(LogLevel messageLevel, const char* message, ...); + static void Log(LogLevel messageLevel, std::string message); + static void Log(std::string message, LogLevel messageLevel = LogLevel::INFO); + static bool IsLoggingEnabled(void); + static LogLevel GetLogLevel(void); + + private: + // Methods + static std::string CreateDateString(void); + static std::string CreateTimestamp(bool appendMS = true, bool iso = true); + static bool CreateDirectory(const std::string& path); + static std::string ConvertLogLevelToString(LogLevel level); + static LogLevel ConvertStringToLogLevel(const std::string& logLevel); + static bool DirectoryExists(const std::string& path); + static std::string GetLogFilePath(void); + static bool LogFileReady(bool appendMode=true); + static void SetLogFilePath(void); + static void SetLoggingFlag(void); + static void SetLogLevel(void); + static void SetLogModuleName(void); + static void SetLogPreferences(void); + static std::string ToUpper(const std::string& input); + static std::string TrimString(const std::string& str); + + // Variables + // Declaring these static so they exist in the class without needing to instantiate it. + static bool loggerInitialized; + static std::string logFilePath; + static bool loggingEnabled; + static std::fstream logFile; + static LogLevel logLevel; + static std::string moduleName; + + static std::shared_ptr loggerInstance; +}; + +#endif diff --git a/include/all.hxx b/include/all.hxx index a3c3fd2..cc64d3b 100755 --- a/include/all.hxx +++ b/include/all.hxx @@ -24,6 +24,7 @@ #include #include #include +#include using namespace std; @@ -53,6 +54,18 @@ struct wetting_front bool to_bottom; // TRUE iff this wetting front is in contact with the layer bottom double dzdt_cm_per_h; // use to store the calculated wetting front speed struct wetting_front *next; // pointer to the next wetting front. + + template + void serialize(Archive &ar, const unsigned int version) { + ar & this->depth_cm; + ar & this->theta; + ar & this->psi_cm; + ar & this->K_cm_per_h; + ar & this->layer_num; + ar & this->front_num; + ar & this->to_bottom; + ar & this->dzdt_cm_per_h; + } }; /* head is a GLOBALLY defined pointer to the first link in the wetting front list. diff --git a/include/bmi_lgar.hxx b/include/bmi_lgar.hxx index b29339a..7c843c5 100644 --- a/include/bmi_lgar.hxx +++ b/include/bmi_lgar.hxx @@ -14,6 +14,8 @@ using namespace std; #include "../bmi/bmi.hxx" #include "all.hxx" #include +#include "vecbuf.hpp" +#include extern "C" { #include "../giuh/giuh.h" @@ -30,7 +32,11 @@ class NotImplemented : public std::logic_error { class BmiLGAR : public bmi::Bmi { public: ~BmiLGAR(); - BmiLGAR():giuh_ordinates(nullptr), giuh_runoff_queue(nullptr) { + BmiLGAR() : + giuh_ordinates(nullptr), + giuh_runoff_queue(nullptr), + m_serialized{} + { this->input_var_names[0] = "precipitation_rate"; this->input_var_names[1] = "potential_evapotranspiration_rate"; this->input_var_names[2] = "soil_temperature_profile"; @@ -130,8 +136,9 @@ public: void global_mass_balance(); double update_calibratable_parameters(); struct model_state* get_model(); - + private: + friend class boost::serialization::access; void realloc_soil(); struct model_state* state; static const int input_var_name_count = 3; @@ -146,6 +153,9 @@ private: double *giuh_ordinates; double *giuh_runoff_queue; + vecbuf m_serialized; + uint64_t m_serialized_length; // needs a stable anchor for GetValuePtr + // unit conversion //struct unit_conversion units; struct bmi_unit_conversion { @@ -163,6 +173,15 @@ private: }; struct bmi_unit_conversion bmi_unit_conv; + + template + void serialize(Archive& ar, const unsigned int version); + template + void serialize_wetting_front_list(Archive &ar, wetting_front **head, int &count); + + void new_serialized(); + void load_serialized(const char* data); + void free_serialized(); }; diff --git a/include/vecbuf.hpp b/include/vecbuf.hpp new file mode 100644 index 0000000..6db982c --- /dev/null +++ b/include/vecbuf.hpp @@ -0,0 +1,115 @@ +#ifndef HPP_STRING_VECBUF +#define HPP_STRING_VECBUF +// https://gist.github.com/stephanlachnit/4a06f8475afd144e73235e2a2584b000 +// SPDX-FileCopyrightText: 2023 Stephan Lachnit +// SPDX-License-Identifier: MIT + +#include +#include +#include + +template> +class vecbuf : public std::basic_streambuf { +public: + using streambuf = std::basic_streambuf; + using char_type = typename streambuf::char_type; + using int_type = typename streambuf::int_type; + using traits_type = typename streambuf::traits_type; + using vector = std::vector; + using value_type = typename vector::value_type; + using size_type = typename vector::size_type; + + // Constructor for vecbuf with optional initial capacity + vecbuf(size_type capacity = 0) : vector_() { reserve(capacity); } + + // Forwarder for std::vector::shrink_to_fit() + constexpr void shrink_to_fit() { vector_.shrink_to_fit(); } + + // Forwarder for std::vector::clear() + constexpr void clear() { vector_.clear(); } + + // Forwarder for std::vector::reserve + constexpr void reserve(size_type capacity) { vector_.reserve(capacity); setp_from_vector(); } + + // Increase the capacity of the buffer by reserving the current_size + additional_capacity + constexpr void reserve_additional(size_type additional_capacity) { reserve(size() + additional_capacity); } + + // Forwarder for std::vector::data + constexpr const value_type* data() const { return vector_.data(); } + + // Forwarder for std::vector::size + constexpr size_type size() const { return vector_.size(); } + + // Forwarder for std::vector::capacity + constexpr size_type capacity() const { return vector_.capacity(); } + + // Implements std::basic_streambuf::xsputn + std::streamsize xsputn(const char_type* s, std::streamsize count) override { + try { + reserve_additional(count); + } + catch (const std::bad_alloc& error) { + // reserve did not work, use slow algorithm + return xsputn_slow(s, count); + } + // reserve worked, use fast algorithm + return xsputn_fast(s, count); + } + +protected: + // Calculates value to std::basic_streambuf::pbase from vector + constexpr value_type* pbase_from_vector() const { return const_cast(vector_.data()); } + + // Calculates value to std::basic_streambuf::pptr from vector + constexpr value_type* pptr_from_vector() const { return const_cast(vector_.data() + vector_.size()); } + + // Calculates value to std::basic_streambuf::epptr from vector + constexpr value_type* epptr_from_vector() const { return const_cast(vector_.data()) + vector_.capacity(); } + + // Sets the values for std::basic_streambuf::pbase, std::basic_streambuf::pptr and std::basic_streambuf::epptr from vector + constexpr void setp_from_vector() { streambuf::setp(pbase_from_vector(), epptr_from_vector()); streambuf::pbump(size()); } + +private: + // std::vector containing the data + vector vector_; + + // Fast implementation of std::basic_streambuf::xsputn if reserve_additional(count) succeeded + std::streamsize xsputn_fast(const char_type* s, std::streamsize count) { + // store current pptr (end of vector location) + auto* old_pptr = pptr_from_vector(); + // resize the vector, does not move since space already reserved + vector_.resize(vector_.size() + count); + // directly memcpy new content to old pptr (end of vector before it was resized) + traits_type::copy(old_pptr, s, count); + // reserve() already calls setp_from_vector(), only adjust pptr to new epptr + streambuf::pbump(count); + + return count; + } + + // Slow implementation of std::basic_streambuf::xsputn if reserve_additional(count) did not succeed, might calls std::basic_streambuf::overflow() + std::streamsize xsputn_slow(const char_type* s, std::streamsize count) { + // reserving entire vector failed, emplace char for char + std::streamsize written = 0; + while (written < count) { + try { + // copy one char, should throw eventually std::bad_alloc + vector_.emplace_back(s[written]); + } + catch (const std::bad_alloc& error) { + // try overflow(), if eof return, else continue + int_type c = this->overflow(traits_type::to_int_type(s[written])); + if (traits_type::eq_int_type(c, traits_type::eof())) { + return written; + } + } + // update pbase, pptr and epptr + setp_from_vector(); + written++; + } + return written; + } + +}; + +#endif diff --git a/src/Logger.cpp b/src/Logger.cpp new file mode 100644 index 0000000..edb7cc3 --- /dev/null +++ b/src/Logger.cpp @@ -0,0 +1,420 @@ +#include "../include/Logger.hpp" +#include +#include +#include +#include // For getenv() +#include +#include // For file handling +#include +#include +#include +#include // For std::string +#include +#include +#include +#include +#include +#include +#include + +const std::string MODULE_NAME = "LASAM"; +const std::string LOG_DIR_NGENCERF = "/ngencerf/data"; // ngenCERF log directory string if environement var empty. +const std::string LOG_DIR_DEFAULT = "run-logs"; // Default parent log directory string if env var empty & ngencerf dosn't exist +const std::string LOG_FILE_EXT = "log"; // Log file name extension +const std::string DS = "/"; // Directory separator +const unsigned int LOG_MODULE_NAME_LEN = 8; // Width of module name for log entries + +const std::string EV_EWTS_LOGGING = "NGEN_EWTS_LOGGING"; // Enable/disable of Error Warning and Trapping System +const std::string EV_NGEN_LOGFILEPATH = "NGEN_LOG_FILE_PATH"; // ngen log file +const std::string EV_MODULE_LOGLEVEL = "LASAM_LOGLEVEL"; // This modules log level +const std::string EV_MODULE_LOGFILEPATH = "LASAM_LOGFILEPATH"; // This modules log full log filename + +bool Logger::loggerInitialized = false; +std::string Logger::logFilePath = ""; +bool Logger::loggingEnabled = true; +std::fstream Logger::logFile; +LogLevel Logger::logLevel = LogLevel::INFO; // Default Log Level +std::string Logger::moduleName = ""; + +std::shared_ptr Logger::loggerInstance; + +using namespace std; + +// String to LogLevel map +static const std::unordered_map logLevelMap = { + {"NONE", LogLevel::NONE }, + {"0", LogLevel::NONE }, + {"DEBUG", LogLevel::DEBUG}, + {"1", LogLevel::DEBUG}, + {"INFO", LogLevel::INFO }, + {"2", LogLevel::INFO }, + {"WARNING", LogLevel::WARNING}, + {"3", LogLevel::WARNING}, + {"SEVERE", LogLevel::SEVERE }, + {"4", LogLevel::SEVERE }, + {"FATAL", LogLevel::FATAL}, + {"5", LogLevel::FATAL}, +}; + +// Reverse map: LogLevel to String +static const std::unordered_map logLevelToStringMap = { + {LogLevel::NONE, "NONE "}, + {LogLevel::DEBUG, "DEBUG "}, + {LogLevel::INFO, "INFO "}, + {LogLevel::WARNING, "WARNING"}, + {LogLevel::SEVERE, "SEVERE "}, + {LogLevel::FATAL, "FATAL "}, +}; + +std::string Logger::ToUpper(const std::string& input) { + std::string result = input; + std::transform(result.begin(), result.end(), result.begin(), + [](unsigned char c){ return std::toupper(c); }); + return result; +} + +// Function to trim leading and trailing spaces +std::string Logger::TrimString(const std::string& str) { + // Trim leading spaces + size_t first = str.find_first_not_of(" \t\n\r\f\v"); + if (first == std::string::npos) { + return ""; // No non-whitespace characters + } + + // Trim trailing spaces + size_t last = str.find_last_not_of(" \t\n\r\f\v"); + + // Return the trimmed string + return str.substr(first, last - first + 1); +} + +LogLevel Logger::ConvertStringToLogLevel(const std::string& levelStr) { + std::string level = TrimString(levelStr); + if (!level.empty()) { + // Convert string to LogLevel (supports both names and numbers) + auto it = logLevelMap.find(level); + if (it != logLevelMap.end()) { + return it->second; // Found valid named or numeric log level + } + + // Try parsing as an integer (for cases where an invalid numeric value is given) + try { + int levelNum = std::stoi(level); + if (levelNum >= 0 && levelNum <= 5) { + return static_cast(levelNum); + } + } catch (...) { + // Ignore errors (e.g., if std::stoi fails for non-numeric input) + } + } + return LogLevel::NONE; +} + +/** + * Convert LogLevel to String Representation of Log Level + * @param logLevel : LogLevel + * @return String log level + */ +std::string Logger::ConvertLogLevelToString(LogLevel level) { + auto it = logLevelToStringMap.find(level); + if (it != logLevelToStringMap.end()) { + return it->second; // Found valid named or numeric log level + } + return "NONE"; +} + +bool Logger::DirectoryExists(const std::string& path) { + struct stat info; + if (stat(path.c_str(), &info) != 0) { + return false; // Cannot access path + } + return (info.st_mode & S_IFDIR) != 0; +} + +/** + * Create the directory checking both the call + * to execute the command and the result of the command + */ +bool Logger::CreateDirectory(const std::string& path) { + + if (!DirectoryExists(path)) { + std::string mkdir_cmd = "mkdir -p " + path; + int status = system(mkdir_cmd.c_str()); + + if (status == -1) { + std::cerr << "system() failed to run mkdir.\n"; + return false; + } else if (WIFEXITED(status)) { + int exitCode = WEXITSTATUS(status); + if (exitCode != 0) { + std::cerr << "mkdir command failed with exit code: " << exitCode << "\n"; + return false; + } + } else { + std::cerr << "mkdir terminated abnormally.\n"; + return false; + } + } + return true; +} + +/** + * Open log file and return open status. If already open, + * ensure the write pointer is at the end of the file. + * + * return bool true if open and good, false otherwise + */ +bool Logger::LogFileReady(bool appendMode) { + + if (logFile.is_open() && logFile.good()) { + logFile.seekp(0, std::ios::end); // Ensure write pointer is at the actual file end + return true; + } else { + // Attempt to open + if (!logFilePath.empty()) { + logFile.open( + logFilePath, + ios::out | ((appendMode) ? ios::app : ios::trunc) + ); // This will silently fail if already open. + if (logFile.good()) { + return true; + } + } + } + return false; +} + +/** + * Set the log file path name using the following pattern + * - Use the module log file if available (unset when first run by ngen), otherwise + * - Use ngen log file if available, otherwise + * - Use /ngencerf/data/run-logs//_ if available, otherwise + * - Use ~/run-logs//_ + * - Onced opened, save the full log path to the modules log environment variable so + * it is only opened once for each ngen run (vs for each catchment) + */ +void Logger::SetLogFilePath(void) { + // Use module logfile path environment variable if it exists + logFilePath = ""; + bool appendEntries = true; + bool mdduleLogEnvExists = false; + const char* envVar = std::getenv(EV_MODULE_LOGFILEPATH.c_str()); // Set once module has successfully opened a log file + if (envVar != nullptr && envVar[0] != '\0') { + logFilePath = envVar; + mdduleLogEnvExists = true; + } else { + envVar = std::getenv(EV_NGEN_LOGFILEPATH.c_str()); // Currently set by ngen-cal but envision set for WCOSS at some point + if (envVar != nullptr && envVar[0] != '\0') { + logFilePath = envVar; + } else { + appendEntries = false; + // Get parent log directory + std::string logFileDir; + if (DirectoryExists(LOG_DIR_NGENCERF)) { + logFileDir = LOG_DIR_NGENCERF + DS + LOG_DIR_DEFAULT; + } else { + logFileDir = "~" + DS + LOG_DIR_DEFAULT; + } + + // Ensure parent log direcotry exists + if (CreateDirectory(logFileDir)) { + // Get full log directory path + const char* envUsername = std::getenv("USER"); + if (envUsername) { + std::string username = envUsername; + logFileDir = logFileDir + DS + username; + } else { + logFileDir = logFileDir + DS + CreateDateString(); + } + // Set the full path if log directory exists/created + if (CreateDirectory(logFileDir)) + logFilePath = logFileDir + DS + MODULE_NAME + "_" + CreateTimestamp(false, false) + + "." + LOG_FILE_EXT; + } + } + } + + // Open log file. If path not found, it will try an alternate file. If neither + // works false returned and logs will be written to stdout. + if (LogFileReady(appendEntries)) { + if (!mdduleLogEnvExists) { + setenv(EV_MODULE_LOGFILEPATH.c_str(), logFilePath.c_str(), 1); + std::cout << "Module " << MODULE_NAME << " Log File: " << logFilePath << std::endl; + LogLevel saveLevel = logLevel; + logLevel = LogLevel::INFO; // Ensure this INFO message is always logged + Log(logLevel, "Logging started. Log File Path: %s\n", logFilePath.c_str()); + logLevel = saveLevel; + } + } else { + std::cout << "Unable to open log file "; + if (!logFilePath.empty()) { + std::cout << logFilePath; + std::cout << " (Perhaps check permissions)" << std::endl; + } + std::cout << "Log entries will be writen to stdout" << std::endl; + } +} + +void Logger::SetLogLevel(void) { + // Set the logger log level if environment var not found + const char* envLogLevel = std::getenv(EV_MODULE_LOGLEVEL.c_str()); + if (envLogLevel != nullptr && envLogLevel[0] != '\0') { + logLevel = ConvertStringToLogLevel(envLogLevel); + } + std::string llMsg = "Log level set to " + ConvertLogLevelToString(logLevel) + "\n"; + cout << MODULE_NAME << " " << llMsg; + LogLevel saveLevel = logLevel; + logLevel = LogLevel::INFO; // Ensure this INFO message is always logged + Log(logLevel, llMsg); + logLevel = saveLevel; + +} + +void Logger::SetLoggingFlag(void) { + const char* envVar = std::getenv(EV_EWTS_LOGGING.c_str()); // Set once module has successfully opened a log file + if (envVar != nullptr && envVar[0] != '\0') { + std::string logState = ToUpper(TrimString(envVar)); + loggingEnabled = (logState == "ENABLED")?true:false; + } + std::cout << MODULE_NAME << " Logging " << ((loggingEnabled)?"ENABLED":"DISABLED") << std::endl; +} + +bool Logger::IsLoggingEnabled(void) { + return loggingEnabled; +} + +void Logger::SetLogModuleName(void) { + // Make sure the module name used for logging entries is all uppercase and 8 characters wide. + moduleName = MODULE_NAME; + std::string upperName = moduleName.substr(0, LOG_MODULE_NAME_LEN); // Truncate to LOG_MODULE_NAME_LEN chars max + std::transform(upperName.begin(), upperName.end(), upperName.begin(), ::toupper); + + std::ostringstream oss; + oss << std::left << std::setw(8) << std::setfill(' ') << upperName; + moduleName = oss.str(); +} + +/** + * Configure Logger Preferences + * @param logFile + * @param level: LogLevel::INFO by Default + * @return void + */ +void Logger::SetLogPreferences(void) { + + if (!loggerInitialized) { + loggerInitialized = true; // Only call this once + + SetLoggingFlag(); + if (loggingEnabled) { + SetLogModuleName(); + SetLogFilePath(); + SetLogLevel(); + } + } +} + +void Logger::Log(LogLevel messageLevel, const char* message, ...) { + va_list args; + va_start(args, message); + + // Make a copy to calculate required size + va_list args_copy; + va_copy(args_copy, args); + int requiredLen = vsnprintf(nullptr, 0, message, args_copy); + va_end(args_copy); + + if (requiredLen > 0) { + std::vector buffer(requiredLen + 1); // +1 for null terminator + vsnprintf(buffer.data(), buffer.size(), message, args); + + va_end(args); + + Log(std::string(buffer.data()), messageLevel); + } else { + va_end(args); // still need to clean up + } +} + +/** + * Log given message with defined parameters and generate message to pass on Console or File + * @param message: Log Message + * @param messageLevel: Log Level, LogLevel::INFO by default + */ +void Logger::Log(LogLevel messageLevel, std::string message) { + Log(message, messageLevel); +} +void Logger::Log(std::string message, LogLevel messageLevel) { + + if (!loggerInitialized) SetLogPreferences(); // Cover case where Log is called before setup done + + // don't log if messageLevel < logLevel + if (loggingEnabled && (messageLevel >= logLevel)) { + std::string logType = ConvertLogLevelToString(messageLevel); + std::string logPrefix = CreateTimestamp() + " " + moduleName + " " + logType; + + // Log message, creating individual entries for a multi-line message + std::istringstream logMsg(message); + std::string line; + if (LogFileReady()) { + while (std::getline(logMsg, line)) { + logFile << logPrefix + " " + line << std::endl; + } + logFile.flush(); + } else { + // Log file not found. Write to stdout. + while (std::getline(logMsg, line)) { + cout << logPrefix + " " + line << std::endl; + } + cout << std::flush; + } + } +} + +std::string Logger::CreateDateString(void) { + std::time_t tt = std::time(0); + std::tm* timeinfo = std::gmtime(&tt); // Use std::localtime(&tt) if you want local time + + char buffer[11]; // Enough for "YYYY-MM-DD" + null terminator + std::strftime(buffer, sizeof(buffer), "%F", timeinfo); // %F == %Y-%m-%d + + std::stringstream ss; + ss << buffer; + + return ss.str(); +} + +std::string Logger::CreateTimestamp(bool appendMS, bool iso) { + using namespace std::chrono; + + // Get current time point + auto now = system_clock::now(); + auto now_time_t = system_clock::to_time_t(now); + + // Get milliseconds + auto ms = duration_cast(now.time_since_epoch()) % 1000; + + // Convert to UTC time + std::tm utc_tm; + gmtime_r(&now_time_t, &utc_tm); + + // Format date/time with strftime + char buffer[32]; + if (iso) { + std::strftime(buffer, sizeof(buffer), "%Y-%m-%dT%H:%M:%S", &utc_tm); + } else { + std::strftime(buffer, sizeof(buffer), "%Y%m%dT%H%M%S", &utc_tm); + } + + if (appendMS) { + // Combine with milliseconds + std::ostringstream oss; + oss << buffer << '.' << std::setw(3) << std::setfill('0') << ms.count(); + return oss.str(); + } + return std::string(buffer); +} + +std::string Logger::GetLogFilePath() { + return logFilePath; +} diff --git a/src/aet.cxx b/src/aet.cxx index 03e4bfa..e37f663 100644 --- a/src/aet.cxx +++ b/src/aet.cxx @@ -2,6 +2,7 @@ #define AET_CXX_INCLUDE #include "../include/all.hxx" +#include "../include/Logger.hpp" //################################################################################ /* authors : Fred Ogden and Ahmad Jan and Peter La Follette @@ -12,7 +13,7 @@ h is the capillary head at the surface and h_50 is the capillary head at which AET = 0.5 * PET. */ //################################################################################ - +char aet_buf[256] = ""; extern double calc_aet(double PET_timestep_cm, double time_step_h, double wilting_point_psi_cm, double field_capacity_psi_cm, int *soil_type, double AET_thresh_Theta, double AET_expon, @@ -20,8 +21,8 @@ extern double calc_aet(double PET_timestep_cm, double time_step_h, double wiltin { if (verbosity.compare("high") == 0) { - printf("Computing AET... \n"); - printf("Note: AET_thresh_theta = %lf and AET_expon = %lf are not used in the computation of the current AET model. \n", AET_thresh_Theta, AET_expon); + LOG(LogLevel::DEBUG,"Computing AET... \n"); + LOG(LogLevel::DEBUG,"Note: AET_thresh_theta = %lf and AET_expon = %lf are not used in the computation of the current AET model. \n", AET_thresh_Theta, AET_expon); } double actual_ET_demand = 0.0; @@ -66,7 +67,7 @@ extern double calc_aet(double PET_timestep_cm, double time_step_h, double wiltin actual_ET_demand = PET_timestep_cm*time_step_h; if (verbosity.compare("high") == 0) { - printf("AET = %14.10f \n",actual_ET_demand); + LOG(LogLevel::DEBUG,"AET = %14.10f \n",actual_ET_demand); } return actual_ET_demand; diff --git a/src/bmi_lgar.cxx b/src/bmi_lgar.cxx index 38e42a1..cee10e3 100644 --- a/src/bmi_lgar.cxx +++ b/src/bmi_lgar.cxx @@ -12,7 +12,12 @@ #include "../bmi/bmi.hxx" #include "../include/bmi_lgar.hxx" #include "../include/all.hxx" +#include "../include/Logger.hpp" +#include +#include + +std::stringstream bmilgar_ss(""); // default verbosity is set to 'none' other option 'high' or 'low' needs to be specified in the config file string verbosity="none"; @@ -32,6 +37,7 @@ BmiLGAR::~BmiLGAR(){ void BmiLGAR:: Initialize (std::string config_file) { + LOG("Inside BmiLGAR::Initialize \n", LogLevel::INFO); if (config_file.compare("") != 0 ) { this->state = new model_state; state->head = NULL; @@ -79,9 +85,10 @@ void BmiLGAR:: Update() { if (verbosity.compare("none") != 0) { - std::cerr<<"---------------------------------------------------------\n"; - std::cerr<<"|****************** LASAM BMI Update... ******************|\n"; - std::cerr<<"---------------------------------------------------------\n"; + bmilgar_ss <<"---------------------------------------------------------\n"; + bmilgar_ss <<"|****************** LASAM BMI Update... ******************|\n"; + bmilgar_ss <<"---------------------------------------------------------\n"; + LOG(bmilgar_ss.str(), LogLevel::INFO); bmilgar_ss.str(""); } double mm_to_cm = 0.1; // unit conversion @@ -180,12 +187,23 @@ Update() double ponded_depth_max_cm = state->lgar_bmi_params.ponded_depth_max_cm; if (verbosity.compare("high") == 0) { - std::cerr<<"Pr [cm/h] (timestep) = "<lgar_bmi_input_params->precipitation_mm_per_h * mm_to_cm <<"\n"; - std::cerr<<"PET [cm/h] (timestep) = "<lgar_bmi_input_params->PET_mm_per_h * mm_to_cm <<"\n"; + bmilgar_ss <<"Pr [cm/h] (timestep) = "<lgar_bmi_input_params->precipitation_mm_per_h * mm_to_cm <<"\n"; + bmilgar_ss <<"PET [cm/h] (timestep) = "<lgar_bmi_input_params->PET_mm_per_h * mm_to_cm <<"\n"; + LOG(bmilgar_ss.str(), LogLevel::INFO); bmilgar_ss.str(""); } - assert (state->lgar_bmi_input_params->precipitation_mm_per_h >= 0.0); - assert(state->lgar_bmi_input_params->PET_mm_per_h >=0.0); + if (state->lgar_bmi_input_params->precipitation_mm_per_h < 0.0) { + std::stringstream error_message; + error_message << "Pr [mm/h] is less than 0: " << state->lgar_bmi_input_params->precipitation_mm_per_h; + LOG(LogLevel::SEVERE, error_message.str()); + throw std::runtime_error(error_message.str()); + } + if (state->lgar_bmi_input_params->PET_mm_per_h < 0.0) { + std::stringstream error_message; + error_message << "PET [mm/h] is less than 0: " << state->lgar_bmi_input_params->PET_mm_per_h; + LOG(LogLevel::SEVERE, error_message.str()); + throw std::runtime_error(error_message.str()); + } // adaptive time step is set if (adaptive_timestep) { @@ -204,7 +222,7 @@ Update() subcycles = state->lgar_bmi_params.forcing_interval; if (verbosity.compare("high") == 0) { - printf("time step size in hours: %lf \n", state->lgar_bmi_params.timestep_h); + LOG(LogLevel::DEBUG,"time step size in hours: %lf \n", state->lgar_bmi_params.timestep_h); } // subcycling loop (loop over model's timestep) @@ -214,8 +232,9 @@ Update() this->state->lgar_bmi_params.timesteps ++; if (verbosity.compare("high") == 0 || verbosity.compare("low") == 0) { - std::cerr<<"BMI Update |---------------------------------------------------------------|\n"; - std::cerr<<"BMI Update |Timesteps = "<< state->lgar_bmi_params.timesteps<<", Time [h] = "<state->lgar_bmi_params.time_s / 3600.<<", Subcycle = "<< cycle <<" of "<state_previous != NULL ){ @@ -251,8 +270,9 @@ Update() //using cerr instead of cout due to some cout buffering issues when running in the ngen framework, cerr doesn't buffer so it prints immediately to the sreeen. if (verbosity.compare("high") == 0 || verbosity.compare("low") == 0) { - std::cerr<<"Pr [cm/h], Pr [cm] (subtimestep), subtimestep [h] = "<lgar_bmi_input_params->precipitation_mm_per_h * mm_to_cm <<", "<< precip_subtimestep_cm <<", "<< subtimestep_h<<" ("<head, state->soil_properties); if (verbosity.compare("high") == 0) { - printf("State before moving creating new WF...\n"); + LOG(LogLevel::DEBUG,"State before moving creating new WF...\n"); listPrint(state->head); } @@ -344,7 +365,7 @@ Update() state->lgar_bmi_params.frozen_factor, &state->head, state->soil_properties); if (verbosity.compare("high") == 0) { - printf("State after moving creating new WF...\n"); + LOG(LogLevel::DEBUG,"State after moving creating new WF...\n"); listPrint(state->head); } @@ -357,7 +378,7 @@ Update() volin_timestep_cm += volin_subtimestep_cm; if (verbosity.compare("high") == 0) { - std::cerr<<"New wetting front created...\n"; + LOG(LogLevel::DEBUG,"New wetting front created...\n"); listPrint(state->head); } } @@ -447,14 +468,14 @@ Update() volQ_gw_timestep_cm += volQ_gw_subtimestep_cm; if (verbosity.compare("high") == 0 || verbosity.compare("low") == 0) { - printf("Printing wetting fronts at this subtimestep... \n"); + LOG(LogLevel::DEBUG,"Printing wetting fronts at this subtimestep... \n"); listPrint(state->head); } bool unexpected_local_error = fabs(local_mb) > 1.0E-4 ? true : false; if (verbosity.compare("high") == 0 || verbosity.compare("low") == 0 || unexpected_local_error) { - printf("\nLocal mass balance at this timestep... \n\ + LOG(LogLevel::DEBUG,"\nLocal mass balance at this timestep... \n\ Error = %14.10f \n\ Initial water = %14.10f \n\ Water added = %14.10f \n\ @@ -468,7 +489,7 @@ Update() volend_subtimestep_cm); if (unexpected_local_error) { - printf("Local mass balance (in this timestep) is %14.10f, larger than expected, needs some debugging...\n ",local_mb); + LOG(LogLevel::DEBUG,"Local mass balance (in this timestep) is %14.10f, larger than expected, needs some debugging...\n ",local_mb); abort(); } @@ -477,7 +498,13 @@ Update() // store local mass balance error to the struct state->lgar_mass_balance.local_mass_balance = local_mb; - assert (state->head->depth_cm > 0.0); // check on negative layer depth --> move this to somewhere else AJ (later) + // check on negative layer depth --> move this to somewhere else AJ (later) + if (state->head->depth_cm <= 0.0) { + std::stringstream error_message; + error_message << "Cycle " << cycle << " has a depth less than or equal to 0: " << state->head->depth_cm; + LOG(LogLevel::SEVERE, error_message.str()); + throw std::runtime_error(error_message.str()); + } bool lasam_standalone = true; #ifdef NGEN @@ -508,15 +535,22 @@ Update() // update thickness/depth and soil moisture of wetting fronts (used for state coupling) struct wetting_front *current = state->head; for (int i=0; ilgar_bmi_params.num_wetting_fronts; i++) { - assert (current != NULL); + if (current == NULL) { + std::stringstream error_message; + error_message << "Wetting front at index " << i << " is null."; + LOG(LogLevel::SEVERE, error_message.str()); + throw std::runtime_error(error_message.str()); + } state->lgar_bmi_params.soil_moisture_wetting_fronts[i] = current->theta; state->lgar_bmi_params.soil_depth_wetting_fronts[i] = current->depth_cm * state->units.cm_to_m; current = current->next; - if (verbosity.compare("high") == 0) - std::cerr<<"Wetting fronts (bmi outputs) (depth in meters, theta)= " + if (verbosity.compare("high") == 0) { + bmilgar_ss <<"Wetting fronts (bmi outputs) (depth in meters, theta)= " <lgar_bmi_params.soil_depth_wetting_fronts[i] <<" "<lgar_bmi_params.soil_moisture_wetting_fronts[i]<<"\n"; - } + LOG(bmilgar_ss.str(), LogLevel::INFO); bmilgar_ss.str(""); + } +} // add to mass balance timestep variables state->lgar_mass_balance.volprecip_timestep_cm = precip_timestep_cm; @@ -564,7 +598,11 @@ Update() void BmiLGAR:: UpdateUntil(double t) { - assert (t > 0.0); + if (t <= 0.0) { + const char *error_message = "Time must be greater than 0."; + LOG(LogLevel::SEVERE, error_message); + throw std::invalid_argument(error_message); + } this->Update(); } @@ -594,17 +632,23 @@ update_calibratable_parameters() layer_num = current->layer_num; soil = state->lgar_bmi_params.layer_soil_type[layer_num]; - assert (current != NULL); + if (current == NULL) { + std::stringstream error_message; + error_message << "Wetting front at index " << i << " is null."; + LOG(LogLevel::SEVERE, error_message.str()); + throw std::invalid_argument(error_message.str()); + } if (verbosity.compare("high") == 0 || verbosity.compare("low") == 0) { - std::cerr<<"----------- Calibratable parameters depending on soil layer (initial values) ----------- \n"; - std::cerr<<"| soil_type = "<< soil <<", layer = "<soil_properties[soil].theta_e = state->lgar_calib_params.theta_e[layer_num-1]; @@ -619,14 +663,15 @@ update_calibratable_parameters() state->soil_properties[soil].theta_e, state->soil_properties[soil].theta_r); if (verbosity.compare("high") == 0 || verbosity.compare("low") == 0) { - std::cerr<<"----------- Calibratable parameters depending on soil layer (updated values) ----------- \n"; - std::cerr<<"| soil_type = "<< soil <<", layer = "<next; @@ -634,18 +679,20 @@ update_calibratable_parameters() //next we update the parameters that apply to the whole model domain and do not depend on soil layer if (verbosity.compare("high") == 0 || verbosity.compare("low") == 0) { - std::cerr<<"----------- Calibratable parameters independent of soil layer (initial values) ----------- \n"; - std::cerr<<"field_capacity_psi = " << state->lgar_bmi_params.field_capacity_psi_cm + bmilgar_ss <<"----------- Calibratable parameters independent of soil layer (initial values) ----------- \n"; + bmilgar_ss <<"field_capacity_psi = " << state->lgar_bmi_params.field_capacity_psi_cm <<", ponded_depth_max = " << state->lgar_bmi_params.ponded_depth_max_cm <<"\n"; + LOG(bmilgar_ss.str(), LogLevel::INFO); bmilgar_ss.str(""); } state->lgar_bmi_params.field_capacity_psi_cm = state->lgar_calib_params.field_capacity_psi; state->lgar_bmi_params.ponded_depth_max_cm = state->lgar_calib_params.ponded_depth_max; if (verbosity.compare("high") == 0 || verbosity.compare("low") == 0) { - std::cerr<<"----------- Calibratable parameters independent of soil layer (updated values) ----------- \n"; - std::cerr<<"field_capacity_psi = " << state->lgar_bmi_params.field_capacity_psi_cm + bmilgar_ss <<"----------- Calibratable parameters independent of soil layer (updated values) ----------- \n"; + bmilgar_ss <<"field_capacity_psi = " << state->lgar_bmi_params.field_capacity_psi_cm <<", ponded_depth_max = " << state->lgar_bmi_params.ponded_depth_max_cm <<"\n"; + LOG(bmilgar_ss.str(), LogLevel::INFO); bmilgar_ss.str(""); } if (verbosity.compare("high") == 0) @@ -653,8 +700,10 @@ update_calibratable_parameters() double volstart_after = lgar_calc_mass_bal(state->lgar_bmi_params.cum_layer_thickness_cm, state->head); - if (verbosity.compare("high") == 0 || verbosity.compare("low") == 0) - std::cerr<<"Mass of water (before and after) = "<< volstart_before<<", "<< volstart_after <<"\n"; + if (verbosity.compare("high") == 0 || verbosity.compare("low") == 0) { + bmilgar_ss <<"Mass of water (before and after) = "<< volstart_before<<", "<< volstart_after <<"\n"; + LOG(bmilgar_ss.str(), LogLevel::INFO); bmilgar_ss.str(""); + } return volstart_after - volstart_before; } @@ -687,36 +736,65 @@ Finalize() delete [] state->lgar_bmi_params.frozen_factor; delete state->lgar_bmi_input_params; delete state; + this->state = NULL; } int BmiLGAR:: GetVarGrid(std::string name) { - if (name.compare("soil_storage_model") == 0 || name.compare("soil_num_wetting_fronts") == 0) // int + if ( + name.compare("soil_storage_model") == 0 + || name.compare("soil_num_wetting_fronts") == 0 + || name.compare("serialization_free") == 0 + ) // int return 0; - else if (name.compare("precipitation_rate") == 0 || name.compare("precipitation") == 0) - return 1; - else if (name.compare("potential_evapotranspiration_rate") == 0 - || name.compare("potential_evapotranspiration") == 0 - || name.compare("actual_evapotranspiration") == 0) // double - return 1; - else if (name.compare("surface_runoff") == 0 || name.compare("giuh_runoff") == 0 - || name.compare("soil_storage") == 0 || name.compare("field_capacity") == 0 || name.compare("ponded_depth_max") == 0)// double - return 1; - else if (name.compare("total_discharge") == 0 || name.compare("infiltration") == 0 - || name.compare("percolation") == 0 || name.compare("groundwater_to_stream_recharge") == 0) // double - return 1; - else if (name.compare("mass_balance") == 0) + else if ( + name.compare("precipitation_rate") == 0 + || name.compare("precipitation") == 0 + || name.compare("potential_evapotranspiration_rate") == 0 + || name.compare("potential_evapotranspiration") == 0 + || name.compare("actual_evapotranspiration") == 0 + || name.compare("surface_runoff") == 0 + || name.compare("giuh_runoff") == 0 + || name.compare("soil_storage") == 0 + || name.compare("field_capacity") == 0 + || name.compare("ponded_depth_max") == 0 + || name.compare("total_discharge") == 0 + || name.compare("infiltration") == 0 + || name.compare("percolation") == 0 + || name.compare("groundwater_to_stream_recharge") == 0 + || name.compare("mass_balance") == 0 + ) // double return 1; - else if (name.compare("soil_depth_layers") == 0 || name.compare("smcmax") == 0 || name.compare("smcmin") == 0 - || name.compare("van_genuchten_m") == 0 || name.compare("van_genuchten_alpha") == 0 || name.compare("van_genuchten_n") == 0 - || name.compare("hydraulic_conductivity") == 0) // array of doubles (fixed length) + else if ( + name.compare("soil_depth_layers") == 0 + || name.compare("smcmax") == 0 + || name.compare("smcmin") == 0 + || name.compare("van_genuchten_m") == 0 + || name.compare("van_genuchten_alpha") == 0 + || name.compare("van_genuchten_n") == 0 + || name.compare("hydraulic_conductivity") == 0 + ) // array of doubles (fixed length) return 2; - else if (name.compare("soil_moisture_wetting_fronts") == 0 || name.compare("soil_depth_wetting_fronts") == 0) // array of doubles (dynamic length) + else if ( + name.compare("soil_moisture_wetting_fronts") == 0 + || name.compare("soil_depth_wetting_fronts") == 0 + ) // array of doubles (dynamic length) return 3; - else if (name.compare("soil_temperature_profile") == 0) // array of doubles (fixed and of the size of soil temperature profile) + else if ( + name.compare("soil_temperature_profile") == 0 + ) // array of doubles (fixed and of the size of soil temperature profile) return 4; + else if ( + name.compare("serialization_state") == 0 + ) // char + return 5; + else if ( + name.compare("serialization_create") == 0 + || name.compare("serialization_size") == 0 + ) // uint64_t + return 6; else return -1; } @@ -731,6 +809,10 @@ GetVarType(std::string name) return "int"; else if (var_grid == 1 || var_grid == 2 || var_grid == 3 || var_grid == 4) return "double"; + else if (var_grid == 5) + return "char"; + else if (var_grid == 6) + return "uint64_t"; else return "none"; } @@ -745,6 +827,10 @@ GetVarItemsize(std::string name) return sizeof(int); else if (var_grid == 1 || var_grid == 2 || var_grid == 3 || var_grid == 4) return sizeof(double); + else if (var_grid == 5) + return sizeof(char); + else if (var_grid == 6) + return sizeof(uint64_t); else return 0; } @@ -858,7 +944,7 @@ GetGridRank(const int grid) int BmiLGAR:: GetGridSize(const int grid) { - if (grid == 0 || grid == 1) + if (grid == 0 || grid == 1 || grid == 6) return 1; else if (grid == 2) // number of layers (fixed) return this->state->lgar_bmi_params.num_layers; @@ -866,6 +952,8 @@ GetGridSize(const int grid) return this->state->lgar_bmi_params.num_wetting_fronts; else if (grid == 4) // number of cells (discretized temperature profile, input from SFT) return this->state->lgar_bmi_params.num_cells_temp; + else if (grid == 5) + return this->m_serialized_length; // serialized state length else return -1; } @@ -937,7 +1025,11 @@ GetValuePtr (std::string name) return (void*)&this->state->lgar_calib_params.ponded_depth_max; else if (name.compare("field_capacity") == 0) return (void*)&this->state->lgar_calib_params.field_capacity_psi; - else { + else if (name.compare("serialization_state") == 0) + return (void*)(this->m_serialized.data()); + else if (name.compare("serialization_size") == 0) { + return (void*)(&this->m_serialized_length); + } else { std::stringstream errMsg; errMsg << "variable "<< name << " does not exist"; throw std::runtime_error(errMsg.str()); @@ -974,6 +1066,17 @@ GetValueAtIndices (std::string name, void *dest, int *inds, int len) void BmiLGAR:: SetValue (std::string name, void *src) { + // state serialization exceptions + if (name == "serialization_state") { + this->load_serialized((char*)src); + return; + } else if (name == "serialization_free") { + this->free_serialized(); + return; + } else if (name.compare("serialization_create") == 0) { + this->new_serialized(); + return; + } void * dest = NULL; dest = this->GetValuePtr(name); @@ -1097,7 +1200,8 @@ void BmiLGAR:: GetGridX(const int grid, double *x) { // this is not needed but printing here to avoid compiler warnings - std::cerr<<"GetGridX: "< +void BmiLGAR:: +serialize(Archive& ar, const unsigned int version) { + model_state* state = this->state; + + // GetValuePtr properties + ar & state->lgar_bmi_input_params->precipitation_mm_per_h; + ar & this->bmi_unit_conv.volprecip_timestep_m; + ar & state->lgar_bmi_input_params->PET_mm_per_h; + ar & this->bmi_unit_conv.volPET_timestep_m; + ar & this->bmi_unit_conv.volAET_timestep_m; + ar & this->bmi_unit_conv.volrunoff_giuh_timestep_m; + ar & this->bmi_unit_conv.volend_timestep_m; + ar & this->bmi_unit_conv.volQ_timestep_m; + ar & this->bmi_unit_conv.volin_timestep_m; + ar & this->bmi_unit_conv.volrech_timestep_m; + ar & this->bmi_unit_conv.volQ_gw_timestep_m; + ar & this->bmi_unit_conv.mass_balance_m; + ar & this->bmi_unit_conv.volrunoff_timestep_m; + ar & state->lgar_bmi_params.num_wetting_fronts; + ar & state->lgar_calib_params.ponded_depth_max; + ar & state->lgar_calib_params.field_capacity_psi; + + // end of update + ar & state->lgar_mass_balance.volprecip_timestep_cm; + ar & state->lgar_mass_balance.volin_timestep_cm; + ar & state->lgar_mass_balance.volon_timestep_cm; + ar & state->lgar_mass_balance.volend_timestep_cm; + ar & state->lgar_mass_balance.volAET_timestep_cm; + ar & state->lgar_mass_balance.volrech_timestep_cm; + ar & state->lgar_mass_balance.volrunoff_timestep_cm; + ar & state->lgar_mass_balance.volQ_timestep_cm; + ar & state->lgar_mass_balance.volQ_gw_timestep_cm; + ar & state->lgar_mass_balance.volPET_timestep_cm; + ar & state->lgar_mass_balance.volrunoff_giuh_timestep_cm; + + ar & state->lgar_mass_balance.volprecip_cm; + ar & state->lgar_mass_balance.volin_cm; + ar & state->lgar_mass_balance.volon_cm; + ar & state->lgar_mass_balance.volend_cm; + ar & state->lgar_mass_balance.volAET_cm; + ar & state->lgar_mass_balance.volrech_cm; + ar & state->lgar_mass_balance.volrunoff_cm; + ar & state->lgar_mass_balance.volQ_cm; + ar & state->lgar_mass_balance.volQ_gw_cm; + ar & state->lgar_mass_balance.volPET_cm; + ar & state->lgar_mass_balance.volrunoff_giuh_cm; + ar & state->lgar_mass_balance.volchange_calib_cm; + + ar & boost::serialization::make_array(state->lgar_bmi_params.cum_layer_thickness_cm, state->lgar_bmi_params.num_layers); + + // giuh state + ar & boost::serialization::make_array(this->giuh_runoff_queue, state->lgar_bmi_params.num_giuh_ordinates); + + // in frozen_factor_hydraulic_conductivity + if (state->lgar_bmi_params.sft_coupled) { + ar & boost::serialization::make_array(state->lgar_bmi_params.frozen_factor, state->lgar_bmi_params.num_layers); + } + + // update_calibratable_parameters + ar & state->lgar_bmi_params.field_capacity_psi_cm; + ar & state->lgar_bmi_params.calib_params_flag; + + // may be set in adapative timesteps + if (state->lgar_bmi_params.adaptive_timestep){ + ar & state->lgar_bmi_params.timestep_h; + } + + // how much time has passed since instantiation + ar & state->lgar_bmi_params.time_s; + ar & state->lgar_bmi_params.timesteps; + + // serialization of arbitrarily-lengthed linked-lists + int num_fronts = state->lgar_bmi_params.num_wetting_fronts; + this->serialize_wetting_front_list(ar, &state->head, state->lgar_bmi_params.num_wetting_fronts); + int num_previous = 0; + this->serialize_wetting_front_list(ar, &state->state_previous, num_previous); + + if (Archive::is_loading::value && num_fronts != state->lgar_bmi_params.num_wetting_fronts) { + // reallocate arrays based on new num_wetting_fronts + this->realloc_soil(); + } + ar & boost::serialization::make_array( + state->lgar_bmi_params.soil_moisture_wetting_fronts, state->lgar_bmi_params.num_wetting_fronts + ); + ar & boost::serialization::make_array( + state->lgar_bmi_params.soil_depth_wetting_fronts, state->lgar_bmi_params.num_wetting_fronts + ); + +} + +template +void BmiLGAR::serialize_wetting_front_list(Archive &ar, wetting_front **head, int &count) { + wetting_front *current; + if (Archive::is_saving::value) { + if (count == 0) // assume recalculation needed if 0 + count = listLength(*head); + ar & count; + current = *head; + while (current != NULL) { + ar & (*current); + current = current->next; + } + } else { // loading + ar & count; + listDelete(*head); + *head = NULL; + wetting_front *prior; + for (int i = 0; i < count; ++i) { + current = new wetting_front(); + ar & (*current); + if (i == 0) { + *head = current; + } else { + prior->next = current; + } + prior = current; + } + } +} + +void BmiLGAR::new_serialized() { + this->m_serialized.clear(); + boost::archive::binary_oarchive archive(this->m_serialized); + try { + archive << (*this); + this->m_serialized_length = this->m_serialized.size(); + } catch (const std::exception &e) { + Logger::Log(LogLevel::SEVERE, "Serializing LASAM encountered an error: %s", e.what()); + this->free_serialized(); + throw; + } +} + +void BmiLGAR::load_serialized(const char* data) { + std::stringstream stream(data); + boost::archive::binary_iarchive archive(stream); + try { + archive >> (*this); + } catch (const std::exception &e) { + Logger::Log(LogLevel::SEVERE, "Deserializing LASAM encountered an error: %s", e.what()); + throw; + } + this->free_serialized(); +} + +void BmiLGAR::free_serialized() { + this->m_serialized.clear(); + this->m_serialized.shrink_to_fit(); + this->m_serialized_length = 0; +} + + #endif diff --git a/src/bmi_main_lgar.cxx b/src/bmi_main_lgar.cxx index a106c8b..1dd07a5 100644 --- a/src/bmi_main_lgar.cxx +++ b/src/bmi_main_lgar.cxx @@ -19,6 +19,7 @@ #include "../bmi/bmi.hxx" #include "../include/all.hxx" #include "../include/bmi_lgar.hxx" +#include "../include/Logger.hpp" // module finds forcing file name in the config file std::string GetForcingFile(std::string config_file); @@ -26,6 +27,7 @@ std::string GetForcingFile(std::string config_file); // module read forcings (precipitation and PET) void ReadForcingData(std::string config_file, std::vector& time, std::vector& precip, std::vector& pet); +std::stringstream bmimain_ss(""); #define SUCCESS 0 @@ -37,9 +39,9 @@ int main(int argc, char *argv[]) bool is_IO_supress = false; // if true no output files will be written if (argc != 2) { - printf("Usage: ./build/xlgar CONFIGURATION_FILE \n"); - printf("Run the LASAM (Lumped Arid/semi-aric Model through its BMI with a configuration file.\n"); - printf("Outputs are written to files `variables_data.csv and layers_data.csv`.\n"); + LOG(LogLevel::DEBUG,"Usage: ./build/xlgar CONFIGURATION_FILE \n"); + LOG(LogLevel::DEBUG,"Run the LASAM (Lumped Arid/semi-aric Model through its BMI with a configuration file.\n"); + LOG(LogLevel::DEBUG,"Outputs are written to files `variables_data.csv and layers_data.csv`.\n"); return SUCCESS; } @@ -88,8 +90,9 @@ int main(int argc, char *argv[]) assert (nsteps <= int(PET.size()) ); // assertion to ensure that nsteps are less or equal than the input data if (verbosity.compare("high") == 0 && !is_IO_supress) { - std::cout<<"Variables are written to file : \'data_variables.csv\' \n"; - std::cout<<"Wetting fronts state is written to file : \'data_layers.csv\' \n"; + bmimain_ss <<"Variables are written to file : \'data_variables.csv\' \n"; + bmimain_ss <<"Wetting fronts state is written to file : \'data_layers.csv\' \n"; + LOG(bmimain_ss.str(), LogLevel::INFO); bmimain_ss.str(""); } FILE *outdata_fptr = NULL; @@ -104,7 +107,7 @@ int main(int argc, char *argv[]) for (int j = 0; j < num_output_var; j++) { fprintf(outdata_fptr,"%s",output_var_names[j].c_str()); if (j == num_output_var-1) - fprintf(outdata_fptr,"\n"); + fprintf(outdata_fptr,"\n"); else fprintf(outdata_fptr,","); } @@ -116,9 +119,10 @@ int main(int argc, char *argv[]) for (int i = 0; i < nsteps; i++) { if (verbosity.compare("none") != 0) { - std::cout<<"===============================================================\n"; - std::cout<<"Real time | "<& time, std::ve std::ifstream fp; fp.open(forcing_file); if (!fp) { - cout<<"file "< #include #include +#include "../include/Logger.hpp" + +std::stringstream lgar_ss(""); using namespace std; @@ -97,7 +100,12 @@ extern void lgar_initialize(string config_file, struct model_state *state) struct wetting_front *current = state->head; for (int i=0; ilgar_bmi_params.num_wetting_fronts; i++) { - assert (current != NULL); + if (current == NULL) { + std::stringstream error_message; + error_message << "Wetting front at index " << i << " is null."; + LOG(LogLevel::SEVERE, error_message.str()); + throw std::runtime_error(error_message.str()); + } soil = state->lgar_bmi_params.layer_soil_type[i+1]; @@ -196,8 +204,8 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) if (param_key == "verbosity") { verbosity = param_value; if (verbosity.compare("none") != 0) { - std::cerr<<"Verbosity is set to \' "<lgar_bmi_params.layer_thickness_cm[layer] = vec[layer-1]; - state->lgar_bmi_params.cum_layer_thickness_cm[layer] = state->lgar_bmi_params.cum_layer_thickness_cm[layer-1] + vec[layer-1]; + state->lgar_bmi_params.cum_layer_thickness_cm[layer] = state->lgar_bmi_params.cum_layer_thickness_cm[layer-1] + vec[layer-1]; } state->lgar_bmi_params.num_layers = vec.size(); @@ -278,11 +287,12 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) is_layer_thickness_set = true; if (verbosity.compare("high") == 0) { - std::cerr<<"Number of layers : "<lgar_bmi_params.num_layers<<"\n"; - for (int i=1; i<=state->lgar_bmi_params.num_layers; i++) - std::cerr<<"Thickness, cum. depth : "<lgar_bmi_params.layer_thickness_cm[i]<<" , " - <lgar_bmi_params.cum_layer_thickness_cm[i]<<"\n"; - std::cerr<<" ***** \n"; + lgar_ss <<"Number of layers : "<lgar_bmi_params.num_layers<<"\n"; + for (int i=1; i<=state->lgar_bmi_params.num_layers; i++) { + lgar_ss <<"Thickness, cum. depth [" << i << "]: "<lgar_bmi_params.layer_thickness_cm[i]<<" , " + <lgar_bmi_params.cum_layer_thickness_cm[i]<<"\n"; + LOG(lgar_ss.str(), LogLevel::INFO); lgar_ss.str(""); + } } continue; @@ -305,15 +315,15 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) giuh_ordinates_temp.resize(vec.size()+1); for (unsigned int i=1; i <= vec.size(); i++) - giuh_ordinates_temp[i] = vec[i-1]; + giuh_ordinates_temp[i] = vec[i-1]; is_giuh_ordinates_set = true; if (verbosity.compare("high") == 0) { - for (int i=1; i <= vec.size(); i++) - std::cerr<<"GIUH ordinates (hourly) : "<lgar_bmi_params.num_cells_temp; i++) - std::cerr<<"Soil z (temperature resolution) : "<lgar_bmi_params.soil_temperature_z[i]<<"\n"; - - std::cerr<<" ***** \n"; + for (int i=0; ilgar_bmi_params.num_cells_temp; i++) { + lgar_ss <<"Soil z (temperature resolution) [" << i << "]: "<lgar_bmi_params.soil_temperature_z[i]<<"\n"; + LOG(lgar_ss.str(), LogLevel::INFO); lgar_ss.str(""); + } } continue; @@ -344,8 +354,8 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) is_initial_psi_set = true; if (verbosity.compare("high") == 0) { - std::cerr<<"Initial Psi : "<lgar_bmi_params.initial_psi_cm<<"\n"; - std::cerr<<" ***** \n"; + lgar_ss <<"Initial Psi : "<lgar_bmi_params.initial_psi_cm<<"\n"; + LOG(lgar_ss.str(), LogLevel::INFO); lgar_ss.str(""); } continue; @@ -360,8 +370,8 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) is_soil_params_file_set = true; if (verbosity.compare("high") == 0) { - std::cerr<<"Soil paramaters file : "<lgar_bmi_params.wilting_point_psi_cm<<"\n"; - std::cerr<<" ***** \n"; + lgar_ss <<"Wilting point Psi [cm] : "<lgar_bmi_params.wilting_point_psi_cm<<"\n"; + LOG(lgar_ss.str(), LogLevel::INFO); lgar_ss.str(""); } continue; @@ -382,8 +392,8 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) is_field_capacity_psi_cm_set = true; if (verbosity.compare("high") == 0) { - std::cerr<<"Field capacity Psi [cm] : "<lgar_bmi_params.field_capacity_psi_cm<<"\n"; - std::cerr<<" ***** \n"; + lgar_ss <<"Field capacity Psi [cm] : "<lgar_bmi_params.field_capacity_psi_cm<<"\n"; + LOG(lgar_ss.str(), LogLevel::INFO); lgar_ss.str(""); } continue; @@ -396,7 +406,8 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) state->lgar_bmi_params.use_closed_form_G = true; } else { - std::cerr<<"Invalid option: use_closed_form_G must be true or false. \n"; + lgar_ss <<"Invalid option: use_closed_form_G must be true or false. \n"; + LOG(lgar_ss.str(), LogLevel::WARNING); lgar_ss.str(""); abort(); } @@ -410,7 +421,8 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) state->lgar_bmi_params.adaptive_timestep = true; } else { - std::cerr<<"Invalid option: adaptive_timestep must be true or false. \n"; + lgar_ss <<"Invalid option: adaptive_timestep must be true or false. \n"; + LOG(lgar_ss.str(), LogLevel::WARNING); lgar_ss.str(""); abort(); } @@ -420,21 +432,26 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) state->lgar_bmi_params.timestep_h = stod(param_value); if (param_unit == "[s]" || param_unit == "[sec]" || param_unit == "") // defalut time unit is seconds - state->lgar_bmi_params.timestep_h /= 3600; // convert to hours + state->lgar_bmi_params.timestep_h /= 3600; // convert to hours else if (param_unit == "[min]" || param_unit == "[minute]") - state->lgar_bmi_params.timestep_h /= 60; // convert to hours + state->lgar_bmi_params.timestep_h /= 60; // convert to hours else if (param_unit == "[h]" || param_unit == "[hr]") - state->lgar_bmi_params.timestep_h /= 1.0; // convert to hours + state->lgar_bmi_params.timestep_h /= 1.0; // convert to hours - assert (state->lgar_bmi_params.timestep_h > 0); + if (state->lgar_bmi_params.timestep_h <= 0) { + std::stringstream error_message; + error_message << "The timestep must be greater than 0. Current value in hours: " << state->lgar_bmi_params.timestep_h; + LOG(LogLevel::SEVERE, error_message.str()); + throw std::runtime_error(error_message.str()); + } is_timestep_set = true; state->lgar_bmi_params.minimum_timestep_h = state->lgar_bmi_params.timestep_h; if (verbosity.compare("high") == 0) { - std::cerr<<"Model timestep [hours,seconds]: "<lgar_bmi_params.timestep_h<<" , " - <lgar_bmi_params.timestep_h*3600<<"\n"; - std::cerr<<" ***** \n"; + lgar_ss <<"Model timestep [hours,seconds]: "<lgar_bmi_params.timestep_h<<" , " + <lgar_bmi_params.timestep_h*3600<<"\n"; + LOG(lgar_ss.str(), LogLevel::INFO); lgar_ss.str(""); } continue; @@ -442,21 +459,26 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) else if (param_key == "endtime") { if (param_unit == "[s]" || param_unit == "[sec]" || param_unit == "") // defalut time unit is seconds - state->lgar_bmi_params.endtime_s = stod(param_value); + state->lgar_bmi_params.endtime_s = stod(param_value); else if (param_unit == "[min]" || param_unit == "[minute]") - state->lgar_bmi_params.endtime_s = stod(param_value) * 60.0; + state->lgar_bmi_params.endtime_s = stod(param_value) * 60.0; else if (param_unit == "[h]" || param_unit == "[hr]") - state->lgar_bmi_params.endtime_s = stod(param_value) * 3600.0; + state->lgar_bmi_params.endtime_s = stod(param_value) * 3600.0; else if (param_unit == "[d]" || param_unit == "[day]") - state->lgar_bmi_params.endtime_s = stod(param_value) * 86400.0; + state->lgar_bmi_params.endtime_s = stod(param_value) * 86400.0; - assert (state->lgar_bmi_params.endtime_s > 0); + if (state->lgar_bmi_params.endtime_s <= 0) { + std::stringstream error_message; + error_message << "The end time must be greater than 0. Current value in seconds: " << state->lgar_bmi_params.endtime_s; + LOG(LogLevel::SEVERE, error_message.str()); + throw std::runtime_error(error_message.str()); + } is_endtime_set = true; if (verbosity.compare("high") == 0) { - std::cerr<<"Endtime [days, hours]: "<< state->lgar_bmi_params.endtime_s/86400.0 <<" , " - << state->lgar_bmi_params.endtime_s/3600.0<<"\n"; - std::cerr<<" ***** \n"; + lgar_ss <<"Endtime [days, hours]: "<< state->lgar_bmi_params.endtime_s/86400.0 <<" , " + << state->lgar_bmi_params.endtime_s/3600.0<<"\n"; + LOG(lgar_ss.str(), LogLevel::INFO); lgar_ss.str(""); } continue; @@ -465,31 +487,37 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) state->lgar_bmi_params.forcing_resolution_h = stod(param_value); if (param_unit == "[s]" || param_unit == "[sec]" || param_unit == "") // defalut time unit is seconds - state->lgar_bmi_params.forcing_resolution_h /= 3600; // convert to hours + state->lgar_bmi_params.forcing_resolution_h /= 3600; // convert to hours else if (param_unit == "[min]" || param_unit == "[minute]") - state->lgar_bmi_params.forcing_resolution_h /= 60; // convert to hours + state->lgar_bmi_params.forcing_resolution_h /= 60; // convert to hours else if (param_unit == "[h]" || param_unit == "[hr]") - state->lgar_bmi_params.forcing_resolution_h /= 1.0; // convert to hours + state->lgar_bmi_params.forcing_resolution_h /= 1.0; // convert to hours - assert (state->lgar_bmi_params.forcing_resolution_h > 0); + if (state->lgar_bmi_params.forcing_resolution_h <= 0) { + std::stringstream error_message; + error_message << "The forcing resolution must be greater than 0. Current value in hours: " << state->lgar_bmi_params.forcing_resolution_h; + LOG(LogLevel::SEVERE, error_message.str()); + throw std::runtime_error(error_message.str()); + } is_forcing_resolution_set = true; if (verbosity.compare("high") == 0) { - std::cerr<<"Forcing resolution [hours]: "<lgar_bmi_params.forcing_resolution_h<<"\n"; - std::cerr<<" ***** \n"; + lgar_ss <<"Forcing resolution [hours]: "<lgar_bmi_params.forcing_resolution_h<<"\n"; + LOG(lgar_ss.str(), LogLevel::INFO); lgar_ss.str(""); } continue; } else if (param_key == "sft_coupled") { if (param_value == "true") { - state->lgar_bmi_params.sft_coupled = 1; + state->lgar_bmi_params.sft_coupled = 1; } else if (param_value == "false") { - state->lgar_bmi_params.sft_coupled = 0; // false + state->lgar_bmi_params.sft_coupled = 0; // false } else { - std::cerr<<"Invalid option: sft_coupled must be true or false. \n"; + lgar_ss <<"Invalid option: sft_coupled must be true or false. \n"; + LOG(lgar_ss.str(), LogLevel::WARNING); lgar_ss.str(""); abort(); } @@ -500,21 +528,22 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) is_ponded_depth_max_cm_set = true; if (verbosity.compare("high") == 0) { - std::cerr<<"Maximum ponded depth [cm] : "<lgar_bmi_params.ponded_depth_max_cm<<"\n"; - std::cerr<<" ***** \n"; + lgar_ss <<"Maximum ponded depth [cm] : "<lgar_bmi_params.ponded_depth_max_cm<<"\n"; + LOG(lgar_ss.str(), LogLevel::INFO); lgar_ss.str(""); } continue; } else if (param_key == "calib_params") { if (param_value == "true") { - state->lgar_bmi_params.calib_params_flag = 1; + state->lgar_bmi_params.calib_params_flag = 1; } else if (param_value == "false") { - state->lgar_bmi_params.calib_params_flag = 0; // false + state->lgar_bmi_params.calib_params_flag = 0; // false } else { - std::cerr<<"Invalid option: calib_params must be true or false. \n"; + lgar_ss <<"Invalid option: calib_params must be true or false. \n"; + LOG(lgar_ss.str(), LogLevel::WARNING); lgar_ss.str(""); abort(); } @@ -526,27 +555,28 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) if (verbosity.compare("high") == 0) { std::string flag = state->lgar_bmi_params.use_closed_form_G == true ? "Yes" : "No"; - std::cerr<<"Using closed_form_G? "<< flag <<"\n"; - std::cerr<<" ***** \n"; + lgar_ss <<"Using closed_form_G? "<< flag <<"\n"; + LOG(lgar_ss.str(), LogLevel::INFO); lgar_ss.str(""); } if (verbosity.compare("high") == 0) { std::string flag = state->lgar_bmi_params.sft_coupled == true ? "Yes" : "No"; - std::cerr<<"Coupled to SoilFreezeThaw? "<< flag <<"\n"; - std::cerr<<" ***** \n"; + lgar_ss <<"Coupled to SoilFreezeThaw? "<< flag <<"\n"; + LOG(lgar_ss.str(), LogLevel::INFO); lgar_ss.str(""); } if(!is_max_valid_soil_types_set) state->lgar_bmi_params.num_soil_types = MAX_NUM_SOIL_TYPES; // maximum number of valid soil types defaults to 15 if (verbosity.compare("high") == 0) { - std::cerr<<"Maximum number of soil types: "<lgar_bmi_params.num_soil_types<<"\n"; - std::cerr<<" ***** \n"; + lgar_ss <<"Maximum number of soil types: "<lgar_bmi_params.num_soil_types<<"\n"; + LOG(lgar_ss.str(), LogLevel::INFO); lgar_ss.str(""); } if (!is_layer_soil_type_set) { stringstream errMsg; errMsg << "The configuration file \'" << config_file <<"\' does not set layer_soil_type. \n"; + LOG(errMsg.str(), LogLevel::WARNING); throw runtime_error(errMsg.str()); } @@ -568,70 +598,79 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) //assert (state->lgar_bmi_params.layer_soil_type[layer] <= state->lgar_bmi_params.num_soil_types); //assert (state->lgar_bmi_params.layer_soil_type[layer] <= max_num_soil_in_file); if (state->lgar_bmi_params.layer_soil_type[layer] > state->lgar_bmi_params.num_soil_types) { - state->lgar_bmi_params.is_invalid_soil_type = true; - if (verbosity.compare("high") == 0) { - std::cerr << "Invalid soil type: " - << state->lgar_bmi_params.layer_soil_type[layer] - <<". Model returns input_precip = ouput_Qout. \n"; - } - break; + state->lgar_bmi_params.is_invalid_soil_type = true; + if (verbosity.compare("high") == 0) { + lgar_ss << "Invalid soil type: " + << state->lgar_bmi_params.layer_soil_type[layer] + <<". Model returns input_precip = ouput_Qout. \n"; + LOG(lgar_ss.str(), LogLevel::WARNING); lgar_ss.str(""); + } + break; } } if (verbosity.compare("high") == 0) { for (int layer=1; layer<=state->lgar_bmi_params.num_layers; layer++) { - int soil = state->lgar_bmi_params.layer_soil_type[layer]; - std::cerr<<"Soil type/name : "<lgar_bmi_params.layer_soil_type[layer] - <<" "<soil_properties[soil].soil_name<<"\n"; + int soil = state->lgar_bmi_params.layer_soil_type[layer]; + lgar_ss <<"Soil type/name : "<lgar_bmi_params.layer_soil_type[layer] + <<" "<soil_properties[soil].soil_name<<"\n"; + LOG(lgar_ss.str(), LogLevel::INFO); lgar_ss.str(""); } - std::cerr<<" ***** \n"; } } else { stringstream errMsg; errMsg << "The configuration file \'" << config_file <<"\' does not set soil_params_file. \n"; + LOG(errMsg.str(), LogLevel::WARNING); throw runtime_error(errMsg.str()); } if (!is_layer_thickness_set) { stringstream errMsg; errMsg << "The configuration file \'" << config_file <<"\' does not set layer_thickness. \n"; + LOG(errMsg.str(), LogLevel::WARNING); throw runtime_error(errMsg.str()); } if (!is_initial_psi_set) { stringstream errMsg; errMsg << "The configuration file \'" << config_file <<"\' does not set initial_psi. \n"; + LOG(errMsg.str(), LogLevel::WARNING); throw runtime_error(errMsg.str()); } if (!is_timestep_set) { stringstream errMsg; errMsg << "The configuration file \'" << config_file <<"\' does not set timestep. \n"; + LOG(errMsg.str(), LogLevel::WARNING); throw runtime_error(errMsg.str()); } if (!is_endtime_set) { stringstream errMsg; errMsg << "The configuration file \'" << config_file <<"\' does not set endtime. \n"; + LOG(errMsg.str(), LogLevel::WARNING); throw runtime_error(errMsg.str()); } if(!is_wilting_point_psi_cm_set) { stringstream errMsg; errMsg << "The configuration file \'" << config_file <<"\' does not set wilting_point_psi. \n Recommended value of 15495.0[cm], corresponding to 15 atm. \n"; + LOG(errMsg.str(), LogLevel::WARNING); throw runtime_error(errMsg.str()); } if(!is_field_capacity_psi_cm_set) { stringstream errMsg; errMsg << "The configuration file \'" << config_file <<"\' does not set field_capacity_psi. \n Recommended value of 340.9[cm] for most soils, corresponding to 1/3 atm, or 103.3[cm] for sands, corresponding to 1/10 atm. \n"; + LOG(errMsg.str(), LogLevel::WARNING); throw runtime_error(errMsg.str()); } if (!is_forcing_resolution_set) { stringstream errMsg; errMsg << "The configuration file \'" << config_file <<"\' does not set forcing_resolution. \n"; + LOG(errMsg.str(), LogLevel::WARNING); throw runtime_error(errMsg.str()); } @@ -649,10 +688,10 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) } if (verbosity.compare("high") == 0) { - for (int i=1; i<=state->lgar_bmi_params.num_giuh_ordinates; i++) - std::cerr<<"GIUH ordinates (scaled) : "<lgar_bmi_params.giuh_ordinates[i]<<"\n"; - - std::cerr<<" ***** \n"; + for (int i=1; i<=state->lgar_bmi_params.num_giuh_ordinates; i++) { + lgar_ss <<"GIUH ordinates (scaled) ["<lgar_bmi_params.giuh_ordinates[i]<<"\n"; + LOG(lgar_ss.str(), LogLevel::INFO); lgar_ss.str(""); + } } giuh_ordinates_temp.clear(); } @@ -660,6 +699,7 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) else if (!is_giuh_ordinates_set) { stringstream errMsg; errMsg << "The configuration file \'" << config_file <<"\' does not set giuh_ordinates. \n"; + LOG(errMsg.str(), LogLevel::WARNING); throw runtime_error(errMsg.str()); } @@ -668,6 +708,7 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) if (!is_soil_z_set) { stringstream errMsg; errMsg << "The configuration file \'" << config_file <<"\' does not set soil_z. \n"; + LOG(errMsg.str(), LogLevel::WARNING); throw runtime_error(errMsg.str()); } } @@ -694,9 +735,9 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) state->lgar_bmi_params.frozen_factor, &state->head, state->soil_properties); if (verbosity.compare("none") != 0) { - std::cerr<<"--- Initial state/conditions --- \n"; + LOG(LogLevel::DEBUG,"--- Initial state/conditions --- \n"); listPrint(state->head); - std::cerr<<" ***** \n"; + LOG(LogLevel::DEBUG," ***** \n"); } // initial mass in the system @@ -706,11 +747,18 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) state->lgar_bmi_params.nint = 120; // hacked, not needed to be an input option state->lgar_bmi_params.num_wetting_fronts = state->lgar_bmi_params.num_layers; - assert (state->lgar_bmi_params.num_layers == listLength(state->head)); + if (state->lgar_bmi_params.num_layers != listLength(state->head)) { + std::stringstream error_message; + error_message << "The number of wetting fronts in the configuration (" << state->lgar_bmi_params.num_layers + << ") is not the same as the number of wetting fronts created in memory (" << listLength(state->head) << ")."; + LOG(LogLevel::SEVERE, error_message.str()); + throw std::runtime_error(error_message.str()); + } if (verbosity.compare("high") == 0) { - std::cerr<<"Initial ponded depth is set to zero. \n"; - std::cerr<<"No. of spatial intervals used in trapezoidal integration to compute G : "<lgar_bmi_params.nint<<"\n"; + lgar_ss <<"Initial ponded depth is set to zero. \n"; + lgar_ss <<"No. of spatial intervals used in trapezoidal integration to compute G : "<lgar_bmi_params.nint<<"\n"; + LOG(lgar_ss.str(), LogLevel::INFO); lgar_ss.str(""); } state->lgar_bmi_input_params = new lgar_bmi_input_parameters; @@ -718,8 +766,9 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) state->lgar_bmi_params.timesteps = 0.0; if (verbosity.compare("none") != 0) { - std::cerr<<"------------- Initialization done! ---------------------- \n"; - std::cerr<<"--------------------------------------------------------- \n"; + lgar_ss <<"------------- Initialization done! ---------------------- \n"; + lgar_ss <<"--------------------------------------------------------- \n"; + LOG(lgar_ss.str(), LogLevel::INFO); lgar_ss.str(""); } } @@ -749,7 +798,7 @@ extern void InitializeWettingFronts(int num_layers, double initial_psi_cm, int * soil_properties[soil].theta_e,soil_properties[soil].theta_r); if (verbosity.compare("high") == 0) { - printf("layer, theta, psi, alpha, m, n, theta_e, theta_r = %d, %6.6f, %6.6f, %6.6f, %6.6f, %6.6f, %6.6f, %6.6f \n", + LOG(LogLevel::DEBUG,"layer, theta, psi, alpha, m, n, theta_e, theta_r = %d, %6.6f, %6.6f, %6.6f, %6.6f, %6.6f, %6.6f, %6.6f \n", layer, theta_init, initial_psi_cm, soil_properties[soil].vg_alpha_per_cm, soil_properties[soil].vg_m, soil_properties[soil].vg_n,soil_properties[soil].theta_e,soil_properties[soil].theta_r); } @@ -833,9 +882,12 @@ extern void frozen_factor_hydraulic_conductivity(struct lgar_bmi_parameters lgar break; } - - // assert (layer_temp > 100.0); // just a check to ensure the while loop executes at least once - assert (count > 0); // just a check to ensure the while loop executes at least once + if (count <= 0) { // just a check to ensure the while loop executes at least once + std::stringstream error_message; + error_message << "The count of layers used for getting layer " << layer << "'s temperature for frozen factor hydraulic conductivity is 0."; + LOG(LogLevel::SEVERE, error_message.str()); + throw std::runtime_error(error_message.str()); + } layer_temp /= count; // layer-averaged temperature @@ -849,8 +901,10 @@ extern void frozen_factor_hydraulic_conductivity(struct lgar_bmi_parameters lgar } if (verbosity.compare("high") == 0) { - for (int i=1; i <= lgar_bmi_params.num_layers; i++) - std::cerr<<"frozen factor = "<< lgar_bmi_params.frozen_factor[i]<<"\n"; + for (int i=1; i <= lgar_bmi_params.num_layers; i++) { + lgar_ss <<"frozen factor ["<0) { @@ -1048,7 +1108,7 @@ extern void lgar_move_wetting_fronts(double timestep_h, double *volin_cm, int wf layer_num_below = (wf == last_wetting_front_index) ? layer_num + 1 : next->layer_num; if (verbosity.compare("high") == 0) { - printf ("Layers (current, above, below) == %d %d %d \n", layer_num, layer_num_above, layer_num_below); + LOG(LogLevel::DEBUG,"Layers (current, above, below) == %d %d %d \n", layer_num, layer_num_above, layer_num_below); listPrint(*head); } @@ -1074,7 +1134,7 @@ extern void lgar_move_wetting_fronts(double timestep_h, double *volin_cm, int wf if ( (wf < last_wetting_front_index) && (layer_num_below != layer_num) ) { if (verbosity.compare("high") == 0) { - printf("case (deepest wetting front within layer) : layer_num (%d) != layer_num_below (%d) \n", layer_num, layer_num_below); + LOG(LogLevel::DEBUG,"case (deepest wetting front within layer) : layer_num (%d) != layer_num_below (%d) \n", layer_num, layer_num_below); } current->theta = calc_theta_from_h(next->psi_cm, vg_a,vg_m, vg_n, theta_e, theta_r); @@ -1098,7 +1158,7 @@ extern void lgar_move_wetting_fronts(double timestep_h, double *volin_cm, int wf if (wf == number_of_wetting_fronts && layer_num_below != layer_num && number_of_wetting_fronts == num_layers) { if (verbosity.compare("high") == 0) { - printf("case (number_of_wetting_fronts equal to num_layers) : l (%d) == num_layers (%d) == num_wetting_fronts(%d) \n", wf, num_layers,number_of_wetting_fronts); + LOG(LogLevel::DEBUG,"case (number_of_wetting_fronts equal to num_layers) : l (%d) == num_layers (%d) == num_wetting_fronts(%d) \n", wf, num_layers,number_of_wetting_fronts); } // local variables @@ -1185,7 +1245,7 @@ extern void lgar_move_wetting_fronts(double timestep_h, double *volin_cm, int wf if (verbosity.compare("high") == 0) { - printf("case (wetting front within a layer) : layer_num (%d) == layer_num_below (%d) \n", layer_num,layer_num_below); + LOG(LogLevel::DEBUG,"case (wetting front within a layer) : layer_num (%d) == layer_num_below (%d) \n", layer_num,layer_num_below); } // if wetting front is the most surficial wetting front @@ -1333,8 +1393,6 @@ extern void lgar_move_wetting_fronts(double timestep_h, double *volin_cm, int wf struct wetting_front *wf_free_drainage = listFindFront(wf_free_drainage_demand, *head, NULL); double mass_timestep = (old_mass + precip_mass_to_add) - (actual_ET_demand + free_drainage_demand); - - assert (old_mass > 0.0); if (fabs(wf_free_drainage->theta - theta_e_k1) < 1E-15) { @@ -1422,7 +1480,7 @@ extern void lgar_move_wetting_fronts(double timestep_h, double *volin_cm, int wf if (verbosity.compare("high") == 0) { - printf("State after moving but before merging wetting fronts...\n"); + LOG(LogLevel::DEBUG,"State after moving but before merging wetting fronts...\n"); listPrint(*head); } @@ -1483,7 +1541,7 @@ extern void lgar_move_wetting_fronts(double timestep_h, double *volin_cm, int wf lgar_fix_dry_over_wet_wetting_fronts(&mass_change, cum_layer_thickness_cm, soil_type, head, soil_properties); if (verbosity.compare("high") == 0) { - printf ("mass change/adjustment (dry_over_wet case) = %lf \n", mass_change); + LOG(LogLevel::DEBUG,"mass change/adjustment (dry_over_wet case) = %lf \n", mass_change); } /* sometimes merging can cause a slight mass balance error, to close the mass balance, the tiny error is @@ -1529,7 +1587,7 @@ extern void lgar_move_wetting_fronts(double timestep_h, double *volin_cm, int wf if (verbosity.compare("high") == 0) - printf("Moving/merging wetting fronts done... \n"); + LOG(LogLevel::DEBUG,"Moving/merging wetting fronts done... \n"); //Just a check to make sure that, when there is only 1 layer, than the existing wetting front is at the correct depth. @@ -1560,9 +1618,9 @@ extern void lgar_merge_wetting_fronts(int *soil_type, double *frozen_factor, str current = *head; if (verbosity.compare("high") == 0) { - printf("State before merging wetting fronts...\n"); + LOG(LogLevel::DEBUG,"State before merging wetting fronts...\n"); listPrint(*head); - printf("Merging wetting fronts... \n"); + LOG(LogLevel::DEBUG,"Merging wetting fronts... \n"); } // local variables @@ -1574,7 +1632,7 @@ extern void lgar_merge_wetting_fronts(int *soil_type, double *frozen_factor, str for (int wf=1; wf != listLength(*head); wf++) { if (verbosity.compare("high") == 0) { - printf("Merge | ********* Wetting Front = %d *********\n", wf); + LOG(LogLevel::DEBUG,"Merge | ********* Wetting Front = %d *********\n", wf); } @@ -1591,7 +1649,12 @@ extern void lgar_merge_wetting_fronts(int *soil_type, double *frozen_factor, str double current_mass_this_layer = current->depth_cm * (current->theta - next->theta) + next->depth_cm*(next->theta - next_to_next->theta); current->depth_cm = current_mass_this_layer / (current->theta - next_to_next->theta); - assert (current->depth_cm > 0.0); + if (current->depth_cm <= 0.0) { + std::stringstream error_message; + error_message << "Wetting front at index " << wf << " has a depth less than or equal to 0. Current value: " << current->depth_cm; + LOG(LogLevel::SEVERE, error_message.str()); + throw std::runtime_error(error_message.str()); + } layer_num = current->layer_num; soil_num = soil_type[layer_num]; @@ -1608,14 +1671,14 @@ extern void lgar_merge_wetting_fronts(int *soil_type, double *frozen_factor, str current->K_cm_per_h = calc_K_from_Se(Se, Ksat_cm_per_h, vg_m); if (verbosity.compare("high") == 0) { - printf ("Deleting wetting front (before)... \n"); + LOG(LogLevel::DEBUG,"Deleting wetting front (before)... \n"); listPrint(*head); } listDeleteFront(next->front_num, head); if (verbosity.compare("high") == 0) { - printf ("Deleting wetting front (after) ... \n"); + LOG(LogLevel::DEBUG,"Deleting wetting front (after) ... \n"); listPrint(*head); } } @@ -1624,7 +1687,7 @@ extern void lgar_merge_wetting_fronts(int *soil_type, double *frozen_factor, str } if (verbosity.compare("high") == 0) { - printf("State after merging wetting fronts...\n"); + LOG(LogLevel::DEBUG,"State after merging wetting fronts...\n"); listPrint(*head); } @@ -1649,13 +1712,13 @@ extern void lgar_wetting_fronts_cross_layer_boundary(int num_layers, current = *head; if (verbosity.compare("high") == 0) { - printf("Layer boundary crossing... \n"); + LOG(LogLevel::DEBUG,"Layer boundary crossing... \n"); } for (int wf=1; wf != listLength(*head); wf++) { if (verbosity.compare("high") == 0) { - printf("Boundary Crossing | ******* Wetting Front = %d ****** \n", wf); + LOG(LogLevel::DEBUG,"Boundary Crossing | ******* Wetting Front = %d ****** \n", wf); } // local variables @@ -1720,7 +1783,7 @@ extern void lgar_wetting_fronts_cross_layer_boundary(int num_layers, } if (verbosity.compare("high") == 0) { - printf("States after wetting fronts cross layer boundary...\n"); + LOG(LogLevel::DEBUG,"States after wetting fronts cross layer boundary...\n"); listPrint(*head); } current = current->next; @@ -1747,7 +1810,7 @@ extern double lgar_wetting_front_cross_domain_boundary(double domain_depth_cm, i int length = listLength(*head); if (verbosity.compare("high") == 0) { - printf("Domain boundary crossing (bottom flux calc.) \n"); + LOG(LogLevel::DEBUG,"Domain boundary crossing (bottom flux calc.) \n"); } // local variables @@ -1759,7 +1822,7 @@ extern double lgar_wetting_front_cross_domain_boundary(double domain_depth_cm, i for (int wf=1; wf != length; wf++) { if (verbosity.compare("high") == 0) { - printf("Domain boundary crossing | ***** Wetting Front = %d ****** \n", wf); + LOG(LogLevel::DEBUG,"Domain boundary crossing | ***** Wetting Front = %d ****** \n", wf); } // ensure that loop iterations never exceed the total number of wetting fronts after altering the list @@ -1803,9 +1866,9 @@ extern double lgar_wetting_front_cross_domain_boundary(double domain_depth_cm, i } if (verbosity.compare("high") == 0) { - printf("State after lowest wetting front contributes to flux through the bottom boundary...\n"); + LOG(LogLevel::DEBUG,"State after lowest wetting front contributes to flux through the bottom boundary...\n"); listPrint(*head); - printf("Bottom boundary flux = %lf \n",bottom_flux_cm); + LOG(LogLevel::DEBUG,"Bottom boundary flux = %lf \n",bottom_flux_cm); } return bottom_flux_cm; @@ -1822,7 +1885,7 @@ extern void lgar_fix_dry_over_wet_wetting_fronts(double *mass_change, double* cu struct wetting_front** head, struct soil_properties_ *soil_properties) { if (verbosity.compare("high") == 0) { - printf("Fix Dry over Wet Wetting Front... \n"); + LOG(LogLevel::DEBUG,"Fix Dry over Wet Wetting Front... \n"); } struct wetting_front *current; @@ -2044,7 +2107,7 @@ extern double lgar_insert_water(bool use_closed_form_G, int nint, double timeste } // if ( (current_mass)/max_storage > 0.99 ){ - // printf("warning: vadose zone is 99 percent full or greater. If you are using the model in an environment with more precipitation than PET, LGAR is not an appropriate model because its lower boundary condition is no flow. \n "); + // LOG(LogLevel::DEBUG,"warning: vadose zone is 99 percent full or greater. If you are using the model in an environment with more precipitation than PET, LGAR is not an appropriate model because its lower boundary condition is no flow. \n "); // } //turning off for now; should really print like once per model run or so, not every time step double ponded_depth_temp = *ponded_depth_cm; @@ -2279,7 +2342,8 @@ extern int lgar_read_vG_param_file(char const* vG_param_file_name, int num_soil_ { if (verbosity.compare("high") == 0) { - std::cerr<<"Reading van Genuchten parameters files...\n"; + lgar_ss <<"Reading van Genuchten parameters files...\n"; + LOG(lgar_ss.str(), LogLevel::INFO); lgar_ss.str(""); } // local vars @@ -2295,7 +2359,7 @@ extern int lgar_read_vG_param_file(char const* vG_param_file_name, int num_soil_ // open the file if((in_vG_params_fptr=fopen(vG_param_file_name,"r"))==NULL) { - printf("Can't open input file named %s. Program stopped.\n",vG_param_file_name); exit(-1); + LOG(LogLevel::DEBUG,"Can't open input file named %s. Program stopped.\n",vG_param_file_name); exit(-1); } fgets(jstr,255,in_vG_params_fptr); // read the header line and ignore @@ -2307,9 +2371,9 @@ extern int lgar_read_vG_param_file(char const* vG_param_file_name, int num_soil_ length=strlen(soil_name); if(length>MAX_SOIL_NAME_CHARS) { - printf("While reading vG soil parameter file: %s, soil name longer than allowed. Increase MAX_SOIL_NAME_CHARS\n", + LOG(LogLevel::DEBUG,"While reading vG soil parameter file: %s, soil name longer than allowed. Increase MAX_SOIL_NAME_CHARS\n", vG_param_file_name); - printf("Program stopped.\n"); + LOG(LogLevel::DEBUG,"Program stopped.\n"); exit(0); } @@ -2330,7 +2394,15 @@ extern int lgar_read_vG_param_file(char const* vG_param_file_name, int num_soil_ soil_properties[soil].bc_lambda = 2.0 / (p - 3.0); soil_properties[soil].bc_psib_cm = (p + 3.0) * (147.8 + 8.1 * p + 0.092 * p * p) / (2.0 * vg_alpha_per_cm * p * (p - 1.0) * (55.6 + 7.4 * p + p * p)); - assert(0.0 < soil_properties[soil].bc_psib_cm); + if (soil_properties[soil].bc_psib_cm <= 0.0) { + std::stringstream error_message; + error_message << "Brooks & Corey bubbling pressure for soil index " + << soil << " (" << soil_properties[soil].soil_name + << ") must be greater than 0. Current value: " + << soil_properties[soil].bc_psib_cm; + LOG(LogLevel::SEVERE, error_message.str()); + throw std::runtime_error(error_message.str()); + } } else { fprintf(stderr, "ERROR: van Genuchten parameter n must be greater than 1\n"); @@ -2363,7 +2435,8 @@ extern void lgar_dzdt_calc(bool use_closed_form_G, int nint, double h_p, int *so double *frozen_factor, struct wetting_front* head, struct soil_properties_ *soil_properties) { if (verbosity.compare("high") == 0) { - std::cerr<<"Calculating dz/dt .... \n"; + lgar_ss <<"Calculating dz/dt .... \n"; + LOG(lgar_ss.str(), LogLevel::INFO); lgar_ss.str(""); } struct wetting_front* current; @@ -2384,6 +2457,7 @@ extern void lgar_dzdt_calc(bool use_closed_form_G, int nint, double h_p, int *so if(head == NULL) { stringstream errMsg; errMsg << "lgar derivative function called for empty list (no wetting front exists) \n"; + LOG(errMsg.str(), LogLevel::WARNING); throw runtime_error(errMsg.str()); } @@ -2400,9 +2474,9 @@ extern void lgar_dzdt_calc(bool use_closed_form_G, int nint, double h_p, int *so K_cm_per_h = current->K_cm_per_h; // K(theta) if (K_cm_per_h < 0) { - printf("K is negative (layer_num, wf_num, K): %d %d %lf \n", layer_num, current->front_num, K_cm_per_h); + LOG(LogLevel::DEBUG,"K is negative (layer_num, wf_num, K): %d %d %lf \n", layer_num, current->front_num, K_cm_per_h); listPrint(head); - printf("Is your n value very close to 1? Very small n values can cause K to become 0. \n"); + LOG(LogLevel::DEBUG,"Is your n value very close to 1? Very small n values can cause K to become 0. \n"); //The parameter n must physically attain a value greater than 1. However, when n is small, and apparently less than 1.02, sometimes n can make K evaluate to 0, for larger values of psi. //So, checking for K_cm_per_h <= 0 has been replaced by checking if K_cm_per_h is negative. K_cm_per_h should never be negative (although perhaps machine precision could make this occur, although we haven't seen it yet), but mathematically can be 0 in some rare cases. abort(); @@ -2445,7 +2519,7 @@ extern void lgar_dzdt_calc(bool use_closed_form_G, int nint, double h_p, int *so } if(theta1 > theta2) { - printf("Calculating dzdt : theta1 > theta2 = (%lf, %lf) ... aborting \n", theta1, theta2); // this should never happen + LOG(LogLevel::DEBUG,"Calculating dzdt : theta1 > theta2 = (%lf, %lf) ... aborting \n", theta1, theta2); // this should never happen exit(0); } @@ -2604,7 +2678,7 @@ extern double lgar_theta_mass_balance(int layer_num, int soil_num, double psi_cm // 1. enough accuracy, 2. the algorithm can't improve the error further, // 3. avoid infinite loop, 4. handles case where theta is very close to theta_r and convergence might be possible but would be extremely slow // 5. handles a corner case when prior mass is tiny (e.g., <1.E-5) - // printf("A1 = %.20f, %.18f %.18f %.18f %.18f \n ",fabs(psi_cm_loc - psi_cm_loc_prev) , psi_cm_loc, psi_cm_loc_prev, factor, delta_mass); + // LOG(LogLevel::DEBUG,"A1 = %.20f, %.18f %.18f %.18f %.18f \n ",fabs(psi_cm_loc - psi_cm_loc_prev) , psi_cm_loc, psi_cm_loc_prev, factor, delta_mass); if (fabs(psi_cm_loc - psi_cm_loc_prev) < 1E-15 && factor < 1E-13) break; diff --git a/src/linked_list.cxx b/src/linked_list.cxx index 0176d41..90483c6 100755 --- a/src/linked_list.cxx +++ b/src/linked_list.cxx @@ -1,4 +1,5 @@ #include "../include/all.hxx" +#include "../include/Logger.hpp" //##################################################################################### /* authors : Fred Ogden and Ahmad Jan @@ -36,7 +37,7 @@ extern void listDelete(struct wetting_front* head) { while (head != NULL) { struct wetting_front *next = head->next; - free( head ); + delete head; head = next; } } @@ -47,15 +48,15 @@ extern void listDelete(struct wetting_front* head) extern void listPrint(struct wetting_front* head) { struct wetting_front *current = head; - printf("\n[ "); + LOG(LogLevel::DEBUG,"\n[ "); //start from the beginning while(current != NULL) { if (current->next == NULL) - printf("(%lf,%6.14f,%d,%d,%d, %e, %lf %6.14f) ] \n",current->depth_cm, current->theta, current->layer_num, + LOG(LogLevel::DEBUG,"(%lf,%6.14f,%d,%d,%d, %e, %lf %6.14f) ] \n",current->depth_cm, current->theta, current->layer_num, current->front_num, current->to_bottom, current->dzdt_cm_per_h, current->K_cm_per_h, current->psi_cm); else - printf("(%lf,%6.14f,%d,%d,%d, %e, %lf %6.14f)\n",current->depth_cm, current->theta, current->layer_num, + LOG(LogLevel::DEBUG,"(%lf,%6.14f,%d,%d,%d, %e, %lf %6.14f)\n",current->depth_cm, current->theta, current->layer_num, current->front_num, current->to_bottom, current->dzdt_cm_per_h, current->K_cm_per_h, current->psi_cm); current = current->next; } @@ -72,7 +73,7 @@ extern struct wetting_front* listCopy(struct wetting_front* current, struct wett } else { - struct wetting_front* wf = (struct wetting_front*)malloc(sizeof(struct wetting_front)); + struct wetting_front* wf = new wetting_front(); wf->depth_cm = current->depth_cm; wf->theta = current->theta; @@ -98,7 +99,7 @@ extern void listInsertFirst(double depth, double theta, int front_num, int layer { //create a link - struct wetting_front *link = (struct wetting_front*) malloc(sizeof(struct wetting_front)); + struct wetting_front *link = new wetting_front(); link->depth_cm = depth; link->theta = theta; @@ -254,7 +255,7 @@ extern struct wetting_front* listDeleteFront(int front_num, struct wetting_front previous = current->next; } - if( current != NULL ) free( current ); + delete current; current = previous; while(previous != NULL) { // decrement all front numbers @@ -282,8 +283,7 @@ extern struct wetting_front* listInsertFront(double depth, double theta, int new if(new_front_num==1) { // create it //create a link - struct wetting_front *link = (struct wetting_front*) malloc(sizeof(struct wetting_front)); - + struct wetting_front *link = new wetting_front(); link->depth_cm = depth; link->theta = theta; link->front_num = new_front_num; @@ -304,7 +304,7 @@ extern struct wetting_front* listInsertFront(double depth, double theta, int new do { if (previous->front_num == new_front_num-1) { // this is where we want to insert it //create a new link - struct wetting_front *link = (struct wetting_front*) malloc(sizeof(struct wetting_front)); + struct wetting_front *link = new wetting_front(); link->depth_cm = depth; link->theta = theta; @@ -352,7 +352,7 @@ extern struct wetting_front* listInsertFrontAtDepth(int num_layers, double *cum_ if(head == NULL) { // list is empty // Kinda weird. Shouldn't call this function to start a list. Create link in the first position - link = (struct wetting_front*) malloc(sizeof(struct wetting_front)); + link = new wetting_front(); link->depth_cm = depth; link->theta = theta; @@ -377,7 +377,7 @@ extern struct wetting_front* listInsertFrontAtDepth(int num_layers, double *cum_ // yep. It's first //create a link and put it at the beginning of the list - link = (struct wetting_front*) malloc(sizeof(struct wetting_front)); + link = new wetting_front(); link->depth_cm = depth; link->theta = theta; @@ -406,7 +406,7 @@ extern struct wetting_front* listInsertFrontAtDepth(int num_layers, double *cum_ // it is in this interval // create a link - link = (struct wetting_front*) malloc(sizeof(struct wetting_front)); + link = new wetting_front(); link->depth_cm = depth; link->theta = theta; @@ -438,7 +438,7 @@ extern struct wetting_front* listInsertFrontAtDepth(int num_layers, double *cum_ current=current->next; while (current != NULL) { // increment all front numbers - //printf("front number: %3d, current->next=%p\n",current->front_num,current->next); + //LOG(LogLevel::DEBUG,"front number: %3d, current->next=%p\n",current->front_num,current->next); current->front_num++; previous = current; current = current->next; @@ -532,7 +532,7 @@ extern void listSortFrontsByDepth(struct wetting_front *head) // probably not needed, left in as an example of how it's done extern void listReverseOrder(struct wetting_front** head_ref) { - printf("In Reverse order.... \n "); + LOG(LogLevel::DEBUG,"In Reverse order.... \n "); struct wetting_front* prev = NULL; struct wetting_front* current = *head_ref; // notice that it is passed down as a pointer to a pointer, that explains struct wetting_front* next; // why they call this "head_ref", it is a pointer referencing the diff --git a/src/main.cxx b/src/main.cxx index c3e51dd..4e71420 100755 --- a/src/main.cxx +++ b/src/main.cxx @@ -1,5 +1,6 @@ #include "../include/all.hxx" // <--- This header file contains all function prototypes and other global definitions. #include +#include "../include/Logger.hpp" // this is a C++ style comment that tells the compiler to ignore the rest of the line. @@ -231,7 +232,7 @@ AET_thresh_Theta = 0.85; // scaled soil moisture (0-1) above which AET=PET AET_expon = 1.0; // exponent that allows curvature of the rising portion of the Budyko curve // check? wilting_point_psi_cm=initial_psi_cm*101325.0/9810.0*100.0; // assumed 15 atm. wilting_point_psi_cm = 154950/10.; // AJ-correction -// printf("wilting::: %lf \n", wilting_point_psi_cm); +//LOG(LogLevel::DEBUG,"wilting::: %lf \n", wilting_point_psi_cm); load_ICs=FALSE; // set to false iff initial conditions are not specified TODO fixeme iff TRUE /*## EXAMPLE ##*/ @@ -256,9 +257,9 @@ soil_type_by_layer[3]=15; // clay //######################################################################################## // option sanity checks. Detect nonsensical or infeasible options or initializations here if(num_soil_types > (MAX_NUM_SOIL_TYPES -1)) - {printf("Too many soil types specified. Increase MAX_NUM_SOIL_TYPES in all.h . Program stopped.\n");exit(0);} + {LOG(LogLevel::DEBUG,"Too many soil types specified. Increase MAX_NUM_SOIL_TYPES in all.h . Program stopped.\n");exit(0);} if(num_soil_layers > (MAX_NUM_SOIL_LAYERS -1)) - {printf("Too many soil types specified. Increase MAX_NUM_SOIL_LAYERS in code. Program stopped.\n");exit(0);} + {LOG(LogLevel::DEBUG,"Too many soil types specified. Increase MAX_NUM_SOIL_LAYERS in code. Program stopped.\n");exit(0);} //########################################################################################## @@ -322,16 +323,16 @@ ymax=cum_layer_thickness_cm[num_soil_layers]; if(debug_flag==TRUE) // this stuff is all demo/test code. { - printf("Initial linked list of wetting fronts:\n"); - //printf("depth,theta,layer,front,to_bottom,dzdt "); + LOG(LogLevel::DEBUG,"Initial linked list of wetting fronts:\n"); + //LOG(LogLevel::DEBUG,"depth,theta,layer,front,to_bottom,dzdt "); //listPrint(); // note in the next line, using the file print function (fprintf), with arg. "stdout" causes write to screen fprintf(stdout,"Initial depth of water stored in soil: %lf cm\n",volstart); fprintf(stdout,"Should be equal to: %lf cm\n",39.1658320808); // a kind of unit test. fprintf(stdout,"difference: %e cm\n",volstart-39.1658320808); - printf("Initial linked list of wetting fronts:\n"); - printf("after calculating dzdt for %d fronts.\n",num_dzdt_calculated); - printf("depth,theta,layer,front,to_bottom,dzdt "); + LOG(LogLevel::DEBUG,"Initial linked list of wetting fronts:\n"); + LOG(LogLevel::DEBUG,"after calculating dzdt for %d fronts.\n",num_dzdt_calculated); + LOG(LogLevel::DEBUG,"depth,theta,layer,front,to_bottom,dzdt "); listPrint(); } @@ -339,9 +340,9 @@ if(illiterate_flag==FALSE) { // open the output file ffor writing if((outl_fptr=fopen("data_layers.out","w"))==NULL) - {printf("Problem opening output file. Program stopped.\n"); exit(0);} + {LOG(LogLevel::DEBUG,"Problem opening output file. Program stopped.\n"); exit(0);} if((outd_fptr=fopen("data_variables.out","w"))==NULL) - {printf("Problem opening output file. Program stopped.\n"); exit(0);} + {LOG(LogLevel::DEBUG,"Problem opening output file. Program stopped.\n"); exit(0);} // add headings to variables forcing file fprintf(outd_fptr,"precip(mm),AET(mm),runoff(mm),ponded_depth(mm),storage(mm),bottom_flux(mm),mass_bal_err(mm)\n"); } @@ -352,7 +353,7 @@ if(illiterate_flag==FALSE) //if((in_forcing_fptr=fopen("Phillipsburg_data1.txt","r"))==NULL) if((in_forcing_fptr=fopen("forcing/Phillipsburg_data2_PET.txt","r"))==NULL) - {printf("Problem opening forcing_data_syn.txt input file. Program stopped.\n"); exit(0);} + {LOG(LogLevel::DEBUG,"Problem opening forcing_data_syn.txt input file. Program stopped.\n"); exit(0);} forcing_interval=(int) (forcing_resolution_s/(time_step_s)+1.0e-08); // add 1.0e-08 to prevent truncation error if(debug_flag) @@ -362,7 +363,7 @@ forcing_interval=(int) (forcing_resolution_s/(time_step_s)+1.0e-08); // add 1.0e if( ((forcing_interval) % (int)(time_step_s+1.0e-08)) != 0) // here % is the remainder operator { - //printf("forcing data interval (s) must be integer divisible by the time step (s). Program stopped.\n"); + //LOG(LogLevel::DEBUG,"forcing data interval (s) must be integer divisible by the time step (s). Program stopped.\n"); //exit(-1); // -1 is arbitrary, and can be targeted. Usually int functions return 0 if successful. } @@ -412,7 +413,7 @@ volstart=lgar_calc_mass_bal(num_soil_layers,cum_layer_thickness_cm); ponded_depth_cm=0.0; // always start simulation with no ponded depth if(debug_flag) - printf("\n ************ ************ TIME STEPPING NOW ****************** \n"); + LOG(LogLevel::DEBUG,"\n ************ ************ TIME STEPPING NOW ****************** \n"); double mm_s_to_cm_hr = 1. * (1.0/10.); // mm per sec to cm per hour conversion double volend_timestep_cm; @@ -431,15 +432,15 @@ for(time_step_num=0;time_step_num0) abort(); ponded_depth_cm += precip_timestep_cm * time_step_h; int wf_free_drainage_demand = wetting_front_free_drainage(); - //printf("wf_free_drainage_demand = %d \n", wf_free_drainage_demand); + //LOG(LogLevel::DEBUG,"wf_free_drainage_demand = %d \n", wf_free_drainage_demand); //bool surficial_wf_created = false; volin_timestep_cm =0.0; @@ -527,7 +528,7 @@ for(time_step_num=0;time_step_num0.0); double mass_source_to_soil_timestep = 0.0; - //printf ("__**___Surficial front create? = %d %d %lf %lf \n ", create_surficial_front, is_top_wf_saturated, head->theta, theta_e); + //LOG(LogLevel::DEBUG,"__**___Surficial front create? = %d %d %lf %lf \n ", create_surficial_front, is_top_wf_saturated, head->theta, theta_e); if(create_surficial_front && !is_top_wf_saturated) // This means that there is no wetting front in the top layer to accept the water, must create one. @@ -537,7 +538,7 @@ for(time_step_num=0;time_step_numtheta; lgar_create_surfacial_front(&ponded_depth_cm, &volin_timestep_cm, dry_depth, theta1, soil_type_by_layer, soil_properties, cum_layer_thickness_cm, nint, time_step_h); @@ -548,7 +549,7 @@ for(time_step_num=0;time_step_num0) //listPrint(); // abort(); @@ -557,7 +558,7 @@ for(time_step_num=0;time_step_num= 0 && precip_this_timestep_cm >0 && precip_previous_timestep_cm >0) { if (ponded_depth_cm > 0 && !create_surficial_front) { - //printf("--- LGAR MAIN **** insert water........ %lf, %d \n", ponded_depth_cm, create_surficial_front); + //LOG(LogLevel::DEBUG,"--- LGAR MAIN **** insert water........ %lf, %d \n", ponded_depth_cm, create_surficial_front); volrunoff_timestep_cm = lgar_insert_water(&ponded_depth_cm, &volin_timestep_cm, precip_timestep_cm, dry_depth, nint, time_step_h, wf_free_drainage_demand, soil_type_by_layer, soil_properties, cum_layer_thickness_cm); //volin_this_timestep = ponded_depth_cm; @@ -565,13 +566,13 @@ for(time_step_num=0;time_step_num1e-7) { - printf("Local mass balance (in this timestep) is %.6e, larger than expected, needs some debugging...\n ",local_mb); + LOG(LogLevel::DEBUG,"Local mass balance (in this timestep) is %.6e, larger than expected, needs some debugging...\n ",local_mb); abort(); } //if (volrech_timestep_cm > 0) @@ -695,32 +696,32 @@ for(time_step_num=0;time_step_numfront_num, temp->layer_num); +// LOG(LogLevel::DEBUG,"\nDeleted value:"); +// LOG(LogLevel::DEBUG,"(%d,%d) ",temp->front_num, temp->layer_num); // } // -//printf("\nList after deleting all items: "); +//LOG(LogLevel::DEBUG,"\nList after deleting all items: "); //listPrint(); //------------------------------------------------- @@ -753,38 +754,38 @@ for(time_step_num=0;time_step_numfront_num,foundLink->layer_num, +// LOG(LogLevel::DEBUG,"Element found: "); +// LOG(LogLevel::DEBUG,"(front: %d in layer: %d at Z=%8.4lf with theta=%8.7lf) ",foundLink->front_num,foundLink->layer_num, // foundLink->depth,foundLink->theta); -// printf("\n"); +// LOG(LogLevel::DEBUG,"\n"); // } //else // { -// printf("Element not found."); +// LOG(LogLevel::DEBUG,"Element not found."); // } //---------------------------------- // EXAMPLE how to delete a front // //listDeleteFront(3); // argument here is wetting front number -//printf("List after deleting wetting front 3: "); +//LOG(LogLevel::DEBUG,"List after deleting wetting front 3: "); //listPrint(); -//printf("\n"); +//LOG(LogLevel::DEBUG,"\n"); //-------------------------------------------------------------------------------------------- // EXAMPLE how to sort a linked list of fronts by depth from surface down to the wetting front // -// printf("\n"); +//LOG(LogLevel::DEBUG,"\n"); // listSortFrontsByDepth(); // -// printf("List after sorting the layer depth: "); +//LOG(LogLevel::DEBUG,"List after sorting the layer depth: "); // listPrint(); //----------------------------------------------------------------- // EXAMPLE how to sort a linked list in terms of decreasing depth // // listReverseOrder(&head); -// printf("\nList after reversing the list: "); +//LOG(LogLevel::DEBUG,"\nList after reversing the list: "); // listPrint(); // EXAMPLE how to delete all fronts @@ -793,11 +794,11 @@ for(time_step_num=0;time_step_numfront_num,temp->layer_num); +// LOG(LogLevel::DEBUG,"\nDeleted value:"); +// LOG(LogLevel::DEBUG,"(%d,%d) ",temp->front_num,temp->layer_num); // } // -// printf("\nList after deleting all items: "); +//LOG(LogLevel::DEBUG,"\nList after deleting all items: "); // listPrint(); //-------------------------------------- @@ -858,7 +859,7 @@ extern double calc_aet(double PET_timestep_cm, double time_step_h, double wiltin //PET_timestep_cm = PET_mm_per_15min*time_step_s/forcing_resolution_s/10.0; // 10 converts from mm to cm - //printf("PET_timestep = %lf %lf \n", PET_timestep_cm, wilting_point_psi_cm); + //LOG(LogLevel::DEBUG,"PET_timestep = %lf %lf \n", PET_timestep_cm, wilting_point_psi_cm); @@ -875,28 +876,28 @@ extern double calc_aet(double PET_timestep_cm, double time_step_h, double wiltin theta = current->theta; double theta_fc = (theta_e-theta_r)*relative_moisture_at_which_PET_equals_AET+theta_r; - //printf("theta_fc = %lf \n", theta_fc); + //LOG(LogLevel::DEBUG,"theta_fc = %lf \n", theta_fc); double wp_head_theta = calc_theta_from_h(wilting_point_psi_cm, vg_a,vg_m, vg_n, theta_e, theta_r); theta_wp = (theta_fc-wp_head_theta)*1/2+wp_head_theta; // theta_50 in python - //printf("theta_50 = %lf %lf %lf \n", theta_wp, wilting_point_psi_cm,wp_head_theta); + //LOG(LogLevel::DEBUG,"theta_50 = %lf %lf %lf \n", theta_wp, wilting_point_psi_cm,wp_head_theta); Se = calc_Se_from_theta(theta_wp,theta_e,theta_r); double psi_wp_cm = calc_h_from_Se(Se, vg_a, vg_m, vg_n); - //printf("psi_50 = %lf %lf \n", psi_wp_cm,current->psi_cm); + //LOG(LogLevel::DEBUG,"psi_50 = %lf %lf \n", psi_wp_cm,current->psi_cm); double h_ratio = 1+pow(current->psi_cm/psi_wp_cm, 3.0); actual_ET_demand = PET_timestep_cm*(1/h_ratio)*time_step_h; - //printf("AET.. = %lf %lf %lf \n", PET_timestep_cm, actual_ET_demand,(1/h_ratio), time_step_h); + //LOG(LogLevel::DEBUG,"AET.. = %lf %lf %lf \n", PET_timestep_cm, actual_ET_demand,(1/h_ratio), time_step_h); if (actual_ET_demand<0) actual_ET_demand=0.0; else if (actual_ET_demand>(PET_timestep_cm*time_step_h)) actual_ET_demand = PET_timestep_cm*time_step_h; - //printf("AET.... = %0.6e %lf %lf \n", actual_ET_demand,(1/h_ratio), time_step_h); + //LOG(LogLevel::DEBUG,"AET.... = %0.6e %lf %lf \n", actual_ET_demand,(1/h_ratio), time_step_h); return actual_ET_demand; /* do @@ -905,12 +906,12 @@ do { wettest_theta_in_root_zone = current->theta; layer = current->layer_num; - if(layer<=0) {printf("problem in calc_aet. Program stopped.\n"); exit(-1);} // should never happen + if(layer<=0) {LOG(LogLevel::DEBUG,"problem in calc_aet. Program stopped.\n"); exit(-1);} // should never happen soil = soil_type[layer]; } current = current->next; } while (current != NULL); - printf("wettest theta = %lf %d \n", wettest_theta_in_root_zone, root_zone_bottom_layer); +LOG(LogLevel::DEBUG,"wettest theta = %lf %d \n", wettest_theta_in_root_zone, root_zone_bottom_layer); // calculate AET from PET and soil moisture Se = calc_Se_from_theta(wettest_theta_in_root_zone, soil_props[soil].theta_e, soil_props[soil].theta_r); diff --git a/src/mem_funcs.cxx b/src/mem_funcs.cxx index b5f1914..373903f 100755 --- a/src/mem_funcs.cxx +++ b/src/mem_funcs.cxx @@ -1,5 +1,6 @@ #include "../include/all.hxx" #include +#include "../include/Logger.hpp" /*#########################################################################*/ /*#########################################################################*/ @@ -30,7 +31,7 @@ int i,frows,fcols; if ((rows==0)||(cols==0)) { - printf("Error: Attempting to allocate array of size 0\n"); + LOG(LogLevel::DEBUG,"Error: Attempting to allocate array of size 0\n"); exit(0); } @@ -66,7 +67,7 @@ int i,frows,fcols; if ((rows==0)||(cols==0)) { - printf("Error: Attempting to allocate array of size 0\n"); + LOG(LogLevel::DEBUG,"Error: Attempting to allocate array of size 0\n"); exit(0); } @@ -101,7 +102,7 @@ void d_alloc(double **var,int size) *var = (double *)malloc(size * sizeof(double)); if (*var == NULL) { - printf("Problem allocating memory for array in d_alloc.\n"); + LOG(LogLevel::DEBUG,"Problem allocating memory for array in d_alloc.\n"); return; } else memset(*var,0,size*sizeof(double)); @@ -118,7 +119,7 @@ void i_alloc(int **var,int size) *var = (int *)malloc(size * sizeof(int)); if (*var == NULL) { - printf("Problem allocating memory in i_alloc\n"); + LOG(LogLevel::DEBUG,"Problem allocating memory in i_alloc\n"); return; } else memset(*var,0,size*sizeof(int)); @@ -135,7 +136,7 @@ void f_alloc(float **var,int size) *var = (float *)malloc(size * sizeof(float)); if (*var == NULL) { - printf("Problem allocating memory for array in f_alloc.\n"); + LOG(LogLevel::DEBUG,"Problem allocating memory for array in f_alloc.\n"); return; } else memset(*var,0,size*sizeof(float)); diff --git a/src/soil_funcs.cxx b/src/soil_funcs.cxx index 7da1adc..d5fd603 100755 --- a/src/soil_funcs.cxx +++ b/src/soil_funcs.cxx @@ -1,4 +1,5 @@ #include "../include/all.hxx" +#include "../include/Logger.hpp" #include "iostream" /*##################################################*/ /*##################################################*/ @@ -60,10 +61,10 @@ extern double calc_Geff(bool use_closed_form_G, double theta1, double theta2, do if (verbosity.compare("high") == 0) { // debug statements to see if calc_Se_from_h function is working properly Se = calc_Se_from_h(h_i,vg_alpha,vg_m,vg_n); - printf("Se_i = %8.6lf, Se_inverse = %8.6lf\n", Se_i, Se); + LOG(LogLevel::DEBUG,"Se_i = %8.6lf, Se_inverse = %8.6lf\n", Se_i, Se); Se = calc_Se_from_h(h_f,vg_alpha,vg_m,vg_n); - printf("Se_f = %8.6lf, Se_inverse = %8.6lf\n", Se_f, Se); + LOG(LogLevel::DEBUG,"Se_f = %8.6lf, Se_inverse = %8.6lf\n", Se_f, Se); } dh = (h_i-h_f)/(double)nint; @@ -114,7 +115,7 @@ extern double calc_Geff(bool use_closed_form_G, double theta1, double theta2, do Geff = fabs(Geff/Ksat); // by convention Geff is a positive quantity if (verbosity.compare("high") == 0){ - printf ("Capillary suction (G) = %8.6lf \n", Geff); + LOG(LogLevel::DEBUG,"Capillary suction (G) = %8.6lf \n", Geff); } return Geff;