diff --git a/CMakeLists.txt b/CMakeLists.txt index cddb89bce..71486868b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -36,7 +36,7 @@ option(ENABLE_MPI "Enable MPI parallel support" OFF) option(ENABLE_OPENMP "Enable OpenMP support" OFF) option(ENABLE_COVERAGE "Enable code coverage output" OFF) option(ENABLE_MEMCHECK "Enable memory checking in tests" OFF) -option(ENABLE_JSON "Enable json configureation file reading" OFF) +option(ENABLE_JSON "Enable json configureation file reading" ON) option(ENABLE_REGRESSION_TESTS "Enable regression tests against the old pre-processed version of micm" ON) option(BUILD_DOCS "Build the documentation" OFF) diff --git a/include/micm/configure/solver_config.hpp b/include/micm/configure/solver_config.hpp new file mode 100644 index 000000000..fd880c387 --- /dev/null +++ b/include/micm/configure/solver_config.hpp @@ -0,0 +1,492 @@ +/* Copyright (C) 2023 National Center for Atmospheric Research, + * + * SPDX-License-Identifier: Apache-2.0 + * + */ + +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace micm +{ + template + class ThrowPolicy + { + class Exception : public std::exception + { + public: + const char* msg_; + + public: + Exception(const char* msg) + : msg_(msg) + { + } + + virtual const char* what() + { + return msg_; + } + }; + + public: + Object OnError(std::string message) + { + throw Exception(message.c_str()); + } + }; + + template + class NoThrowPolicy + { + public: + Object OnError(std::string message) + { + std::cerr << message << std::endl; + return Object(); + } + }; + + // Solver parameters + struct SolverParameters + { + micm::System system_; + std::vector processes_; + }; + + // Error code + enum class ConfigErrorCode + { + None = 0, + NoAcess, + FileNotFound, + KeyNotFound, + }; + + // JSON Configure paser + template + class JsonReaderPolicy : public ErrorPolicy + { + using json = nlohmann::json; + + public: + std::vector species_; + std::vector emissions_; + std::vector first_order_loss_; + std::vector processes_; + micm::Phase gas_phase_; + std::unordered_map phases_; + + // Constants + static const inline std::string SPECIES_CONFIG = + "species.json"; // TODO:jiwon 6/6 - instead of searching, pass the configure path + static const inline std::string REACTIONS_CONFIG = "mechanism.json"; // TODO:jiwon 6/6 + + static const inline std::string CAMP_FILES = "camp-files"; + static const inline std::string CAMP_DATA = "camp-data"; + + static const inline std::string TYPE = "type"; + + // Functions + + /// @brief read and parse JSON objects + /// @param + /// @return SolverParameters if parsing is success, else returns ConfigErrorCode + std::variant ReadAndParse(const std::filesystem::path& path) + { + // Check whether file exists + if (!std::filesystem::exists(path)) + { + std::string err_msg = "Configuration file at path " + path.string() + " does not exist\n"; + this->OnError(err_msg); + + return micm::ConfigErrorCode::FileNotFound; + } + + // Read file to get the list of configure files + json data = json::parse(std::ifstream(path)); + if (!ValidateJsonWithKey(data, CAMP_FILES)) + { + return micm::ConfigErrorCode::KeyNotFound; + } + + // Check whether the listed files exist and determine the sequence to read files. + std::string species_file; + std::vector other_files; + bool found_species_file = false; + + for (const auto& file : data[CAMP_FILES].get>()) + { + if (!std::filesystem::exists(file)) + { + std::string err_msg = "Configuration file at path " + file + " does not exist\n"; + this->OnError(err_msg); + + return micm::ConfigErrorCode::FileNotFound; + } + + // Find species file to read first + std::size_t found = file.find(SPECIES_CONFIG); + if (found != std::string::npos) + { + species_file = file; + found_species_file = true; + } + else + { + other_files.push_back(file); + } + } + + // Read species file to create Species and Phase that needs to be known to System and Process + if (found_species_file) + { + if (!ConfigureSpecies(species_file)) + { + return micm::ConfigErrorCode::KeyNotFound; + } + } + else + { + std::string err_msg = "Species configure file does not exist\n"; + this->OnError(err_msg); + + return micm::ConfigErrorCode::FileNotFound; + } + + // Read files, eg. reactions.json + for (const auto& file : other_files) + { + json file_data = json::parse(std::ifstream(file)); + if (!ValidateJsonWithKey(file_data, CAMP_DATA)) + { + return micm::ConfigErrorCode::KeyNotFound; + } + + std::vector objects; + for (const auto& element : file_data[CAMP_DATA]) + { + objects.push_back(element); + } + + if (!ParseObjectArray(objects)) + { + return micm::ConfigErrorCode::KeyNotFound; + } + } + + micm::SystemParameters sysParams = { gas_phase_, phases_ }; + + return micm::SolverParameters{ micm::System(sysParams), processes_ }; + } + + /// @brief Create 'Species' and 'Phase' + /// @param path to 'Species' file + /// @return True at success + bool ConfigureSpecies(const std::string& file) + { + json file_data = json::parse(std::ifstream(file)); + + if (!ValidateJsonWithKey(file_data, CAMP_DATA)) + return false; + + std::vector objects; + for (const auto& element : file_data[CAMP_DATA]) + objects.push_back(element); + + for (const auto& object : objects) + { + if (!ValidateJsonWithKey(object, TYPE)) + return false; + + std::string type = object[TYPE].get(); + + if (type == "CHEM_SPEC") + { + if (!ParseChemicalSpecies(object)) + return false; + } + else if (type == "RELATIVE_TOLERANCE") + { + if (!ParseRelativeTolerance(object)) + return false; + } + } + // After creating Species, create Phase + gas_phase_.species_ = species_; + + return true; + } + + bool ValidateJsonWithKey(const json& object, const std::string& key) + { + if (!object.contains(key)) + { + this->OnError("Key " + key + " was not found in the config file"); + + return false; + } + return true; + } + + bool ParseObjectArray(const std::vector& objects) + { + for (const auto& object : objects) + { + if (!ValidateJsonWithKey(object, TYPE)) + return false; + + std::string type = object[TYPE].get(); + + if (type == "MECHANISM") + { + if (!ParseMechanism(object)) + return false; + } + else if (type == "PHOTOLYSIS") + { + if (!ParsePhotolysis(object)) + return false; + } + else if (type == "ARRHENIUS") + { + if (!ParseArrhenius(object)) + return false; + } + else if (type == "EMISSION") + { + if (!ParseEmission(object)) + return false; + } + else if (type == "FIRST_ORDER_LOSS") + { + if (!ParseFirstOrderLoss(object)) + return false; + } + else + { + this->OnError("Unknown key in config file: " + type); + return false; + } + } + return true; + } + + bool ParseChemicalSpecies(const json& object) + { + std::vector required_keys = { "name" }; + std::vector optional_keys = { "absolute tolerance" }; + + for (const auto& key : required_keys) + { + if (!ValidateJsonWithKey(object, key)) + return false; + } + + std::string name = object["name"].get(); + + std::string key = "absolute tolerance"; + + if (object.contains(key)) + { + double abs_tol = object[key].get(); + auto species = Species(name, Property(key, "", abs_tol)); + species_.push_back(species); + } + else + { + species_.push_back(Species(name)); + } + + return true; + } + + bool ParseRelativeTolerance(const json& object) + { + // TODO: what is this? + return true; + } + + bool ParseMechanism(const json& object) + { + std::vector required_keys = { "name", "reactions" }; + for (const auto& key : required_keys) + { + if (!ValidateJsonWithKey(object, key)) + return false; + } + std::vector objects; + for (const auto& element : object["reactions"]) + { + objects.push_back(element); + } + + return ParseObjectArray(objects); + } + + bool ParsePhotolysis(const json& object) + { + const std::string REACTANTS = "reactants"; + const std::string PRODUCTS = "products"; + const std::string MUSICA_NAME = "MUSICA name"; + const std::string YIELD = "yield"; + + const double DEFAULT_YEILD = 1.0; + + for (const auto& key : { REACTANTS, PRODUCTS, MUSICA_NAME }) + { + if (!ValidateJsonWithKey(object, key)) + return false; + } + + // Create process + std::vector reactants; + for (auto& [key, value] : object[REACTANTS].items()) + { + reactants.push_back(micm::Species(key)); + } + + std::vector> products; + for (auto& [key, value] : object[PRODUCTS].items()) + { + if (value.contains(YIELD)) + { + products.push_back(std::make_pair(micm::Species(key), value[YIELD])); + } + else + { + products.push_back(std::make_pair(micm::Species(key), DEFAULT_YEILD)); + } + } + + std::unique_ptr rate_ptr = + std::make_unique(object[MUSICA_NAME].get()); + + processes_.push_back(micm::Process(reactants, products, std::move(rate_ptr), gas_phase_)); + + return true; + } + + bool ParseArrhenius(const json& object) + { + const std::string REACTANTS = "reactants"; + const std::string PRODUCTS = "products"; + const std::string YIELD = "yield"; + + const double DEFAULT_YEILD = 1.0; + + // Check required json objects exist + for (const auto& key : { REACTANTS, PRODUCTS }) + { + if (!ValidateJsonWithKey(object, key)) + return false; + } + + // Create process + std::vector reactants; + for (auto& [key, value] : object[REACTANTS].items()) + { + reactants.push_back(micm::Species(key)); + } + + std::vector> products; + for (auto& [key, value] : object[PRODUCTS].items()) + { + if (value.contains(YIELD)) + { + products.push_back(std::make_pair(micm::Species(key), value[YIELD])); + } + else + { + products.push_back(std::make_pair(micm::Species(key), DEFAULT_YEILD)); + } + } + + micm::ArrheniusRateConstantParameters parameters; + if (object.contains("A")) + { + parameters.A_ = object["A"].get(); + } + if (object.contains("B")) + { + parameters.B_ = object["B"].get(); + } + if (object.contains("C")) + { + parameters.C_ = object["C"].get(); + } + if (object.contains("D")) + { + parameters.D_ = object["D"].get(); + } + if (object.contains("E")) + { + parameters.E_ = object["E"].get(); + } + + std::unique_ptr rate_ptr = + std::make_unique(micm::ArrheniusRateConstantParameters(parameters)); + + processes_.push_back(micm::Process(reactants, products, std::move(rate_ptr), gas_phase_)); + + return true; + } + + bool ParseEmission(const json& object) + { + std::vector required_keys = { "species" }; + std::vector optional_keys = { "MUSICA name" }; + for (const auto& key : required_keys) + { + if (!ValidateJsonWithKey(object, key)) + return false; + } + + std::string name = object["species"].get(); + + emissions_.push_back(Species(name)); + + return true; + } + + bool ParseFirstOrderLoss(const json& object) + { + std::vector required_keys = { "species" }; + std::vector optional_keys = { "MUSICA name" }; + for (const auto& key : required_keys) + { + if (!ValidateJsonWithKey(object, key)) + return false; + } + + std::string name = object["species"].get(); + + first_order_loss_.push_back(Species(name)); + + return true; + } + }; + + template class ConfigTypePolicy = JsonReaderPolicy, template class ErrorPolicy = NoThrowPolicy> + class SolverConfig : public ConfigTypePolicy>> + { + public: + std::variant Configure(const std::filesystem::path& path) + { + return this->ReadAndParse(path); + } + }; + +} // namespace micm \ No newline at end of file diff --git a/include/micm/configure/system_builder.hpp b/include/micm/configure/system_builder.hpp deleted file mode 100644 index 33c2d7b8c..000000000 --- a/include/micm/configure/system_builder.hpp +++ /dev/null @@ -1,319 +0,0 @@ -/* Copyright (C) 2023 National Center for Atmospheric Research, - * - * SPDX-License-Identifier: Apache-2.0 - * - */ - -#pragma once - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#ifdef USE_JSON -# include -#endif - -namespace micm -{ - - template - class ThrowPolicy - { - class Exception : public std::exception - { - public: - const char* msg_; - - public: - Exception(const char* msg) - : msg_(msg) - { - } - - virtual const char* what() - { - return msg_; - } - }; - - public: - Object OnError(std::string message) - { - throw Exception(message.c_str()); - } - }; - - template - class NoThrowPolicy - { - public: - Object OnError(std::string message) - { - std::cerr << message << std::endl; - return Object(); - } - }; - - template - class JsonReaderPolicy : public ErrorPolicy - { - using json = nlohmann::json; - std::vector species_; - std::vector emissions_; - std::vector first_order_loss_; - std::vector> photolysis_reactions_; - std::vector> arrhenius_reactions_; - - public: - std::unique_ptr ReadAndParse(std::filesystem::path path) - { - species_.clear(); - emissions_.clear(); - first_order_loss_.clear(); - photolysis_reactions_.clear(); - arrhenius_reactions_.clear(); - - if (!std::filesystem::exists(path)) - { - std::string err_msg = "Configuration file at path " + path.string() + " does not exist\n"; - return this->OnError(err_msg); - } - - json data = json::parse(std::ifstream(path)); - - std::string key = "camp-files"; - ValidateJsonWithKey(data, key); - - std::vector files = data[key].get>(); - - key = "camp-data"; - for (const auto& file : files) - { - if (!std::filesystem::exists(path)) - { - std::string err_msg = "Configuration file at path " + path.string() + " does not exist\n"; - return this->OnError(err_msg); - } - json file_data = json::parse(std::ifstream(file)); - ValidateJsonWithKey(file_data, key); - std::vector objects; - for (const auto& element : file_data[key]) - { - objects.push_back(element); - } - - ParseObjectArray(objects); - } - - micm::SystemParameters parameters; - parameters.arrhenius_reactions_ = arrhenius_reactions_; - parameters.photolysis_reactions_ = photolysis_reactions_; - parameters.gas_phase_ = micm::Phase(species_); - - return std::make_unique(micm::System(parameters)); - } - - void ValidateJsonWithKey(const json& object, std::string key) - { - if (!object.contains(key)) - { - this->OnError("Key " + key + " was not found in the config file"); - } - } - - void ParseObjectArray(const std::vector& objects) - { - for (const auto& object : objects) - { - std::string key = "type"; - ValidateJsonWithKey(object, key); - std::string type = object[key].get(); - - if (type == "CHEM_SPEC") - { - ParseChemicalSpecies(object); - } - else if (type == "RELATIVE_TOLERANCE") - { - ParseRelativeTolerance(object); - } - else if (type == "MECHANISM") - { - ParseMechanism(object); - } - else if (type == "PHOTOLYSIS") - { - ParsePhotolysis(object); - } - else if (type == "ARRHENIUS") - { - ParseArrhenius(object); - } - else if (type == "EMISSION") - { - ParseEmission(object); - } - else if (type == "FIRST_ORDER_LOSS") - { - ParseFirstOrderLoss(object); - } - else - { - this->OnError("Unknown key in config file: " + type); - } - } - } - - void ParseChemicalSpecies(const json& object) - { - std::vector required_keys = { "name" }; - std::vector optional_keys = { "absolute tolerance" }; - for (const auto& key : required_keys) - ValidateJsonWithKey(object, key); - - std::string name = object["name"].get(); - - std::string key = "absolute tolerance"; - - if (object.contains(key)) - { - double abs_tol = object[key].get(); - auto species = Species(name, Property(key, "", abs_tol)); - species_.push_back(species); - } - else - { - species_.push_back(Species(name)); - } - } - - void ParseRelativeTolerance(const json& object) - { - // TODO: what is this? - } - - void ParseMechanism(const json& object) - { - std::vector required_keys = { "name", "reactions" }; - for (const auto& key : required_keys) - ValidateJsonWithKey(object, key); - - std::vector objects; - for (const auto& element : object["reactions"]) - { - objects.push_back(element); - } - - ParseObjectArray(objects); - } - - void ParsePhotolysis(const json& object) - { - std::vector required_keys = { "reactants", "products", "MUSICA name" }; - for (const auto& key : required_keys) - ValidateJsonWithKey(object, key); - - std::vector reactants; - for (auto& [key, value] : object["reactants"].items()) - { - reactants.push_back(Species(key)); - } - std::vector products; - for (auto& [key, value] : object["products"].items()) - { - products.push_back(Species(key)); - } - std::string name = object["MUSICA name"].get(); - - using reaction = IntraPhaseProcess; - photolysis_reactions_.push_back(reaction(reactants, products, PhotolysisRateConstant(0, name))); - } - - void ParseArrhenius(const json& object) - { - std::vector required_keys = { "reactants", "products" }; - for (const auto& key : required_keys) - ValidateJsonWithKey(object, key); - - std::vector reactants; - for (auto& [key, value] : object["reactants"].items()) - { - reactants.push_back(Species(key)); - } - std::vector products; - for (auto& [key, value] : object["products"].items()) - { - products.push_back(Species(key)); - } - - ArrheniusRateConstantParameters parameters; - - if (object.contains("A")) - { - parameters.A_ = object["A"].get(); - } - if (object.contains("B")) - { - parameters.B_ = object["B"].get(); - } - if (object.contains("C")) - { - parameters.C_ = object["C"].get(); - } - if (object.contains("D")) - { - parameters.D_ = object["D"].get(); - } - if (object.contains("E")) - { - parameters.E_ = object["E"].get(); - } - - using reaction = IntraPhaseProcess; - arrhenius_reactions_.push_back(reaction(reactants, products, ArrheniusRateConstant(parameters))); - } - - void ParseEmission(const json& object) - { - std::vector required_keys = { "species" }; - std::vector optional_keys = { "MUSICA name" }; - for (const auto& key : required_keys) - ValidateJsonWithKey(object, key); - - std::string name = object["species"].get(); - - emissions_.push_back(Species(name)); - } - - void ParseFirstOrderLoss(const json& object) - { - std::vector required_keys = { "species" }; - std::vector optional_keys = { "MUSICA name" }; - for (const auto& key : required_keys) - ValidateJsonWithKey(object, key); - - std::string name = object["species"].get(); - - first_order_loss_.push_back(Species(name)); - } - }; - - template class ConfigTypePolicy = JsonReaderPolicy, template class ErrorPolicy = ThrowPolicy> - class SystemBuilder : public ConfigTypePolicy>> - { - public: - std::unique_ptr Build(std::filesystem::path path) - { - return this->ReadAndParse(path); - } - }; - -} // namespace micm \ No newline at end of file diff --git a/include/micm/process/arrhenius_rate_constant.hpp b/include/micm/process/arrhenius_rate_constant.hpp index e00f5fd26..ed8131a04 100644 --- a/include/micm/process/arrhenius_rate_constant.hpp +++ b/include/micm/process/arrhenius_rate_constant.hpp @@ -9,9 +9,6 @@ namespace micm { - - class System; - struct ArrheniusRateConstantParameters { /// @brief Pre-exponential factor, (cm−3)^(−(𝑛−1))s−1 @@ -38,13 +35,13 @@ namespace micm const ArrheniusRateConstantParameters parameters_; public: - /// @brief Default constructor. All terms will be zero + /// @brief Default constructor ArrheniusRateConstant(); /// @brief An explicit constructor where each term can be set. Set B and E to zero to get the common form of the /// Arrhenius equation /// @param parameters A set of arrhenius rate constants - ArrheniusRateConstant(ArrheniusRateConstantParameters parameters); + ArrheniusRateConstant(const ArrheniusRateConstantParameters& parameters); /// @brief Deep copy std::unique_ptr clone() const override; @@ -53,7 +50,8 @@ namespace micm /// @param conditions The current environmental conditions of the chemical system /// @param custom_parameters User-defined rate constant parameters /// @return A rate constant based off of the conditions in the system - double calculate(const Conditions& conditions, std::vector::const_iterator custom_parameters) const override; + double calculate(const Conditions& conditions, const std::vector::const_iterator& custom_parameters) + const override; double calculate(const double& temperature, const double& pressure) const; }; @@ -63,7 +61,7 @@ namespace micm { } - inline ArrheniusRateConstant::ArrheniusRateConstant(ArrheniusRateConstantParameters parameters) + inline ArrheniusRateConstant::ArrheniusRateConstant(const ArrheniusRateConstantParameters& parameters) : parameters_(parameters) { } @@ -73,8 +71,9 @@ namespace micm return std::unique_ptr{ new ArrheniusRateConstant{ *this } }; } - inline double ArrheniusRateConstant::calculate(const Conditions& conditions, std::vector::const_iterator custom_parameters) - const + inline double ArrheniusRateConstant::calculate( + const Conditions& conditions, + const std::vector::const_iterator& custom_parameters) const { return calculate(conditions.temperature_, conditions.pressure_); } diff --git a/include/micm/process/external_rate_constant.hpp b/include/micm/process/external_rate_constant.hpp deleted file mode 100644 index dd3681618..000000000 --- a/include/micm/process/external_rate_constant.hpp +++ /dev/null @@ -1,46 +0,0 @@ -/* Copyright (C) 2023 National Center for Atmospheric Research, - * - * SPDX-License-Identifier: Apache-2.0 - */ -#pragma once - -#include -#include - -namespace micm -{ - - /** - * @brief A rate constant from an external model - * - */ - class ExternalRateConstant : public RateConstant - { - public: - /// @brief The rate - const double rate_; - /// @brief The condition this rate applies to - const Condition condition_; - - public: - /// @brief Default constructor - ExternalRateConstant(); - - /// @brief Calculate a reaction rate - /// @param system A defined system - /// @return A reaction rate - double calculate(const System& system) override; - }; - - inline ExternalRateConstant::ExternalRateConstant() - : rate_(), - condition_("", "") - { - } - - inline double ExternalRateConstant::calculate(const System& system) - { - return 0.0; - } - -} // namespace micm \ No newline at end of file diff --git a/include/micm/process/intraphase_process.hpp b/include/micm/process/intraphase_process.hpp deleted file mode 100644 index 081a4693d..000000000 --- a/include/micm/process/intraphase_process.hpp +++ /dev/null @@ -1,41 +0,0 @@ -/* Copyright (C) 2023 National Center for Atmospheric Research, - * - * SPDX-License-Identifier: Apache-2.0 - */ -#pragma once - -#include -#include - -namespace micm -{ - - /** - * @brief An intraphase process - * - */ - template - class IntraPhaseProcess - { - public: - const std::vector reactants_; - const std::vector products_; - const Rate rate_; - - public: - /// @brief - /// @param reactants - /// @param products - /// @param rate - IntraPhaseProcess(std::vector reactants, std::vector products, Rate rate); - }; - - template - inline IntraPhaseProcess::IntraPhaseProcess(std::vector reactants, std::vector products, Rate rate) - : reactants_(reactants), - products_(products), - rate_(rate) - { - } - -} // namespace micm \ No newline at end of file diff --git a/include/micm/process/intraphase_process_builder.hpp b/include/micm/process/intraphase_process_builder.hpp deleted file mode 100644 index d4b354a85..000000000 --- a/include/micm/process/intraphase_process_builder.hpp +++ /dev/null @@ -1,86 +0,0 @@ -/* Copyright (C) 2023 National Center for Atmospheric Research, - * - * SPDX-License-Identifier: Apache-2.0 - */ -#pragma once - -#include -#include -#include - -namespace micm -{ - - /** - * @brief Creates any type of intraphase process builder - * - */ - class IntraPhaseProcessBuilder - { - private: - public: - /// @brief Adds an additional phase to the state of the process - /// @param phase A phase - /// @return A reference to this object - IntraPhaseProcessBuilder& For(const Phase& phase); - /// @brief Adds an additional species to the process - /// @param species A species - /// @return A reference to this object - IntraPhaseProcessBuilder& With(const Species& species); - /// @brief Adds a reacting species - /// @param reactant A Species - /// @return A reference to this object - IntraPhaseProcessBuilder& Reacting(const Species& reactant); - /// @brief Adds a species that will be produced - /// @param product A species - /// @return A reference to this object - IntraPhaseProcessBuilder& Producing(const Species& product); - /// @brief Provides a yield in the amount of TODO UNITS - /// @param yield A value - /// @return A reference to this object - IntraPhaseProcessBuilder& WithYield(const double yield); - /// @brief Add a rate contant - /// @param rate_constant A rate constant - /// @return A reference to this object - IntraPhaseProcessBuilder& WithRateConstant(const RateConstant& rate_constant); - /// @brief Create the final process - /// @return A reference to this object - IntraPhaseProcessBuilder& Build(); - }; - - inline IntraPhaseProcessBuilder& IntraPhaseProcessBuilder::For(const Phase& phase) - { - return *this; - } - - inline IntraPhaseProcessBuilder& IntraPhaseProcessBuilder::With(const Species& species) - { - return *this; - } - - inline IntraPhaseProcessBuilder& IntraPhaseProcessBuilder::Reacting(const Species& reactant) - { - return *this; - } - - inline IntraPhaseProcessBuilder& IntraPhaseProcessBuilder::Producing(const Species& product) - { - return *this; - } - - inline IntraPhaseProcessBuilder& IntraPhaseProcessBuilder::WithYield(const double yield) - { - return *this; - } - - inline IntraPhaseProcessBuilder& IntraPhaseProcessBuilder::WithRateConstant(const RateConstant& rate_constant) - { - return *this; - } - - inline IntraPhaseProcessBuilder& IntraPhaseProcessBuilder::Build() - { - return *this; - } - -} // namespace micm \ No newline at end of file diff --git a/include/micm/process/photolysis_rate_constant.hpp b/include/micm/process/photolysis_rate_constant.hpp index 0d36df790..9a64a1a37 100644 --- a/include/micm/process/photolysis_rate_constant.hpp +++ b/include/micm/process/photolysis_rate_constant.hpp @@ -17,7 +17,7 @@ namespace micm class PhotolysisRateConstant : public RateConstant { public: - const std::string name_; + std::string name_; public: /// @brief Default constructor. @@ -25,7 +25,7 @@ namespace micm /// @brief /// @param name A name for this reaction - PhotolysisRateConstant(const std::string name); + PhotolysisRateConstant(const std::string& name); /// @brief Deep copy std::unique_ptr clone() const override; @@ -42,7 +42,8 @@ namespace micm /// @param conditions The current environmental conditions of the chemical system /// @param custom_parameters User-defined rate constant parameters /// @return A rate constant based off of the conditions in the system - double calculate(const Conditions& conditions, std::vector::const_iterator custom_parameters) const override; + double calculate(const Conditions& conditions, const std::vector::const_iterator& custom_parameters) + const override; }; inline PhotolysisRateConstant::PhotolysisRateConstant() @@ -50,7 +51,7 @@ namespace micm { } - inline PhotolysisRateConstant::PhotolysisRateConstant(const std::string name) + inline PhotolysisRateConstant::PhotolysisRateConstant(const std::string& name) : name_(name) { } @@ -60,8 +61,9 @@ namespace micm return std::unique_ptr{ new PhotolysisRateConstant{ *this } }; } - inline double PhotolysisRateConstant::calculate(const Conditions& conditions, std::vector::const_iterator custom_parameters) - const + inline double PhotolysisRateConstant::calculate( + const Conditions& conditions, + const std::vector::const_iterator& custom_parameters) const { return (double)*custom_parameters; } diff --git a/include/micm/process/process.hpp b/include/micm/process/process.hpp index 0cf2b9df9..ce2fcc13a 100644 --- a/include/micm/process/process.hpp +++ b/include/micm/process/process.hpp @@ -22,13 +22,13 @@ namespace micm }; class ProcessBuilder; - + struct Process { - const std::vector reactants_; - const std::vector products_; - const std::unique_ptr rate_constant_; - const Phase phase_; + std::vector reactants_; + std::vector products_; + std::unique_ptr rate_constant_; + Phase phase_; /// @brief Update the solver state rate constants /// @param processes The set of processes being solved @@ -39,6 +39,14 @@ namespace micm static ProcessBuilder create(); Process(ProcessBuilder& builder); Process(const Process& other); + + Process(const std::vector& reactants, const std::vector& products, std::unique_ptr rate_constant, const Phase& phase) + : reactants_(reactants), + products_(products), + rate_constant_(std::move(rate_constant)), + phase_(phase) + { + } }; class ProcessBuilder diff --git a/include/micm/process/rate_constant.hpp b/include/micm/process/rate_constant.hpp index 1d3bf69ae..1243250f2 100644 --- a/include/micm/process/rate_constant.hpp +++ b/include/micm/process/rate_constant.hpp @@ -11,9 +11,6 @@ namespace micm { - - class System; - /** * @brief A base class for any type of rate constant * @@ -27,22 +24,20 @@ namespace micm virtual std::unique_ptr clone() const = 0; /// @brief Returns the number of doubles needed to hold user-defined rate constant parameters /// @return Number of user-defined rate constant parameters - virtual std::size_t SizeCustomParameters() const; + virtual std::size_t SizeCustomParameters() const + { + return 0; + } + /// @brief Calculate the rate constant for a set of conditions /// @param conditions The current environmental conditions of the chemical system /// @param custom_parameters User defined rate constant parameters /// @return The reaction rate constant - virtual double calculate(const Conditions& conditions, std::vector::const_iterator custom_parameters) const + virtual double calculate(const Conditions& conditions, const std::vector::const_iterator& custom_parameters) + const { return 0; - }; - - private: + } }; - inline std::size_t RateConstant::SizeCustomParameters() const - { - return 0; - } - } // namespace micm \ No newline at end of file diff --git a/include/micm/process/troe_rate_constant.hpp b/include/micm/process/troe_rate_constant.hpp index fba14dbb9..cf1319889 100644 --- a/include/micm/process/troe_rate_constant.hpp +++ b/include/micm/process/troe_rate_constant.hpp @@ -40,21 +40,22 @@ namespace micm const TroeRateConstantParameters parameters_; public: - /// @brief Default constructor is not allowed - TroeRateConstant() = delete; + /// @brief Default constructor + TroeRateConstant(); /// @brief An explicit constructor /// @param parameters A set of troe rate constants - TroeRateConstant(TroeRateConstantParameters parameters); + TroeRateConstant(const TroeRateConstantParameters& parameters); + + /// @brief Deep copy + std::unique_ptr clone() const override; /// @brief Calculate the rate constant /// @param conditions The current environmental conditions of the chemical system /// @param custom_parameters User-defined rate constant parameters /// @return A rate constant based off of the conditions in the system - double calculate(const Conditions& conditions, std::vector::const_iterator custom_parameters) const override; - - /// @brief Deep copy - std::unique_ptr clone() const override; + double calculate(const Conditions& conditions, const std::vector::const_iterator& custom_parameters) + const override; /// @brief Calculate the rate constant /// @param temperature Temperature in [K] @@ -63,7 +64,12 @@ namespace micm double calculate(const double& temperature, const double& air_number_density) const; }; - inline TroeRateConstant::TroeRateConstant(TroeRateConstantParameters parameters) + inline TroeRateConstant::TroeRateConstant() + : parameters_() + { + } + + inline TroeRateConstant::TroeRateConstant(const TroeRateConstantParameters& parameters) : parameters_(parameters) { } @@ -73,7 +79,9 @@ namespace micm return std::unique_ptr{ new TroeRateConstant{ *this } }; } - inline double TroeRateConstant::calculate(const Conditions& conditions, std::vector::const_iterator custom_parameters) const + inline double TroeRateConstant::calculate( + const Conditions& conditions, + const std::vector::const_iterator& custom_parameters) const { return calculate(conditions.temperature_, conditions.air_density_); } diff --git a/include/micm/solver/state.hpp b/include/micm/solver/state.hpp index 08e76823b..48af8bd3f 100644 --- a/include/micm/solver/state.hpp +++ b/include/micm/solver/state.hpp @@ -4,6 +4,7 @@ #include #include #include +#include namespace micm { diff --git a/include/micm/system/phase.hpp b/include/micm/system/phase.hpp index 6a79adadf..e54636777 100644 --- a/include/micm/system/phase.hpp +++ b/include/micm/system/phase.hpp @@ -19,24 +19,24 @@ namespace micm { public: /// @brief The list of species - const std::vector species_; + std::vector species_; public: - Phase(); - Phase operator=(const Phase& other); + /// @brief Default constructor + Phase() + : species_() + { + } /// @brief Create a phase with a set of species /// @param species A unique list of species - Phase(const std::vector species) + Phase(const std::vector& species) : species_(species) { } - }; - Phase::Phase() - : species_() - { - } + Phase operator=(const Phase& other); + }; Phase Phase::operator=(const Phase& other) { diff --git a/include/micm/util/sparse_matrix.hpp b/include/micm/util/sparse_matrix.hpp index 5db9e26d8..80bf4e6b3 100644 --- a/include/micm/util/sparse_matrix.hpp +++ b/include/micm/util/sparse_matrix.hpp @@ -97,6 +97,16 @@ namespace micm { } + SparseMatrix& operator=(SparseMatrixBuilder& builder) + { + number_of_blocks_ = builder.number_of_blocks_; + data_ = std::vector(builder.NumberOfElements(), builder.initial_value_); + row_ids_ = builder.RowIdsVector(); + row_start_ = builder.RowStartVector(); + + return *this; + } + std::vector& AsVector() { return data_; diff --git a/test/integration/chapman.cpp b/test/integration/chapman.cpp index 9034900a7..bd742456f 100644 --- a/test/integration/chapman.cpp +++ b/test/integration/chapman.cpp @@ -12,6 +12,42 @@ using yields = std::pair; +#ifdef USE_JSON +# include +TEST(ChapmanIntegration, CanBuildChapmanSystemUsingConfig) +{ + micm::SolverConfig solverConfig{}; // Throw policy + std::variant configs = + solverConfig.Configure("./unit_configs/chapman/config.json"); + + // Check if parsing is successful and returns 'Solverparameters' + auto* solver_params_ptr = std::get_if(&configs); + EXPECT_TRUE(solver_params_ptr != nullptr); + + micm::SolverParameters& solver_params = *solver_params_ptr; + + micm::RosenbrockSolver solver{ solver_params.system_, + std::move(solver_params.processes_), + micm::RosenbrockSolverParameters{} }; + + micm::State state = solver.GetState(); + + std::vector concentrations{ 0.1, 0.1, 0.1, 0.2, 0.2, 0.2, 0.3, 0.3, 0.3 }; + state.variables_[0] = concentrations; + std::vector photo_rates{ 0.1, 0.2, 0.3 }; + state.custom_rate_parameters_[0] = photo_rates; + state.conditions_[0].temperature_ = 2; + state.conditions_[0].pressure_ = 3; + + for (double t{}; t < 100; ++t) + { + state.custom_rate_parameters_[0] = photo_rates; + auto result = solver.Solve(t, t + 0.5, state); + // output state + } +} +#endif + TEST(ChapmanIntegration, CanBuildChapmanSystem) { auto o = micm::Species("O"); @@ -26,33 +62,32 @@ TEST(ChapmanIntegration, CanBuildChapmanSystem) micm::Phase gas_phase{ std::vector{ o, o1d, o2, o3, m, ar, n2, h2o, co2 } }; - micm::Process r1 = - micm::Process::create() - .reactants({ o1d, n2 }) - .products({ yields(o, 1), yields(n2, 1) }) - .rate_constant(micm::ArrheniusRateConstant({ .A_ = 2.15e-11, .B_= 0, .C_ = 110 })) - .phase(gas_phase); - - micm::Process r2 = - micm::Process::create() - .reactants({ o1d, o2 }) - .products({ yields(o, 1), yields(o2, 1) }) - .rate_constant(micm::ArrheniusRateConstant(micm::ArrheniusRateConstantParameters{ .A_ = 3.3e-11, .B_ = 0, .C_ = 55 })) - .phase(gas_phase); - - micm::Process r3 = - micm::Process::create() - .reactants({ o, o3 }) - .products({ yields(o2, 2) }) - .rate_constant(micm::ArrheniusRateConstant(micm::ArrheniusRateConstantParameters{ .A_ = 8e-12, .B_ = 0, .C_ = -2060 })) - .phase(gas_phase); - - micm::Process r4 = - micm::Process::create() - .reactants({ o, o2, m }) - .products({ yields(o3, 1), yields(m, 1) }) - .rate_constant(micm::ArrheniusRateConstant(micm::ArrheniusRateConstantParameters{ .A_ = 6.0e-34, .B_ = 0, .C_ = 2.4 })) - .phase(gas_phase); + micm::Process r1 = micm::Process::create() + .reactants({ o1d, n2 }) + .products({ yields(o, 1), yields(n2, 1) }) + .rate_constant(micm::ArrheniusRateConstant({ .A_ = 2.15e-11, .B_ = 0, .C_ = 110 })) + .phase(gas_phase); + + micm::Process r2 = micm::Process::create() + .reactants({ o1d, o2 }) + .products({ yields(o, 1), yields(o2, 1) }) + .rate_constant(micm::ArrheniusRateConstant( + micm::ArrheniusRateConstantParameters{ .A_ = 3.3e-11, .B_ = 0, .C_ = 55 })) + .phase(gas_phase); + + micm::Process r3 = micm::Process::create() + .reactants({ o, o3 }) + .products({ yields(o2, 2) }) + .rate_constant(micm::ArrheniusRateConstant( + micm::ArrheniusRateConstantParameters{ .A_ = 8e-12, .B_ = 0, .C_ = -2060 })) + .phase(gas_phase); + + micm::Process r4 = micm::Process::create() + .reactants({ o, o2, m }) + .products({ yields(o3, 1), yields(m, 1) }) + .rate_constant(micm::ArrheniusRateConstant( + micm::ArrheniusRateConstantParameters{ .A_ = 6.0e-34, .B_ = 0, .C_ = 2.4 })) + .phase(gas_phase); micm::Process photo_1 = micm::Process::create() .reactants({ o2 }) diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt index b844bc2f4..2f61327bd 100644 --- a/test/unit/CMakeLists.txt +++ b/test/unit/CMakeLists.txt @@ -1,4 +1,6 @@ -add_subdirectory(configure) +if(ENABLE_JSON) + add_subdirectory(configure) +endif() add_subdirectory(process) add_subdirectory(solver) add_subdirectory(system) diff --git a/test/unit/configure/CMakeLists.txt b/test/unit/configure/CMakeLists.txt index 86fc68e93..f8b91a056 100644 --- a/test/unit/configure/CMakeLists.txt +++ b/test/unit/configure/CMakeLists.txt @@ -6,4 +6,4 @@ include(test_util) ################################################################################ # Tests -#create_standard_test(NAME system_builder SOURCES test_system_builder.cpp) +create_standard_test(NAME solver_config SOURCES test_solver_config.cpp) \ No newline at end of file diff --git a/test/unit/configure/test_solver_config.cpp b/test/unit/configure/test_solver_config.cpp new file mode 100644 index 000000000..25d9254d2 --- /dev/null +++ b/test/unit/configure/test_solver_config.cpp @@ -0,0 +1,130 @@ +#include + +#ifdef USE_JSON +# include + +TEST(SolverConfig, DetectsInvalidConfigFileAndThrow) +{ + micm::SolverConfig solverConfig{}; + EXPECT_ANY_THROW(solverConfig.Configure("not_a_config_file.json")); +} + +TEST(SolverConfig, DetectsInvalidConfigFileAndNoThrowDoesntThrow) +{ + micm::SolverConfig solverConfig{}; + EXPECT_NO_THROW(solverConfig.Configure("not_a_config_file.json")); + + std::variant configs = solverConfig.Configure("not_a_config_file.json"); + EXPECT_EQ(std::get(configs), micm::ConfigErrorCode::FileNotFound); +} + +TEST(SolverConfig, ReadAndParseSystem) +{ + micm::SolverConfig solverConfig{}; + std::variant configs = + solverConfig.Configure("./unit_configs/chapman/config.json"); + + // Check if parsing is successful and returns 'Solverparameters' + auto* solver_params_ptr = std::get_if(&configs); + EXPECT_TRUE(solver_params_ptr != nullptr); + + micm::SolverParameters& solver_params = *solver_params_ptr; + + // Check 'gas_phase' in 'System' + EXPECT_EQ(solver_params.system_.gas_phase_.species_.size(), 9); + + // Check 'phases' in 'System' + EXPECT_EQ(solver_params.system_.phases_.size(), 0); + + // Check 'name' and 'properties' in 'Species' + std::vector> species_name_and_num_properties = { + std::make_pair("M", 0), std::make_pair("Ar", 1), std::make_pair("CO2", 1), + std::make_pair("H2O", 1), std::make_pair("N2", 1), std::make_pair("O1D", 1), + std::make_pair("O", 1), std::make_pair("O2", 1), std::make_pair("O3", 1) + }; + + short idx = 0; + for (const auto& s : solver_params.system_.gas_phase_.species_) + { + EXPECT_EQ(s.name_, species_name_and_num_properties[idx].first); + EXPECT_EQ(s.properties_.size(), species_name_and_num_properties[idx].second); + idx++; + } +} + +TEST(SolverConfig, ReadAndParseProcess) +{ + micm::SolverConfig solverConfig{}; + std::variant configs = + solverConfig.Configure("./unit_configs/chapman/config.json"); + + // Check if parsing is successful and returns 'Solverparameters' + auto* solver_params_ptr = std::get_if(&configs); + EXPECT_TRUE(solver_params_ptr != nullptr); + + micm::SolverParameters& solver_params = *solver_params_ptr; + auto& process_vector = solver_params.processes_; + + // Check the number of 'Process' created + EXPECT_EQ(process_vector.size(), 7); + + // Check the number of 'reactants' and 'products' in each 'Process' + // Check 'yield' value for the first product and the number of 'spieces in 'phase' in each 'Process' + int num_reactants_in_each_process[] = { 1, 1, 1, 2, 2, 2, 3 }; + int num_products_in_each_process[] = { 1, 2, 2, 2, 2, 1, 2 }; + double yield_value_of_first_product_in_each_process[] = { 2.0, 1.0, 1.0, 1.0, 1.0, 2.0, 1.0 }; + int num_phase_in_each_process = 9; + + short idx = 0; + for (const auto& p : process_vector) + { + EXPECT_EQ(p.reactants_.size(), num_reactants_in_each_process[idx]); + EXPECT_EQ(p.products_.size(), num_products_in_each_process[idx]); + EXPECT_EQ(p.products_[0].second, yield_value_of_first_product_in_each_process[idx]); + EXPECT_EQ(p.phase_.species_.size(), num_phase_in_each_process); + idx++; + } + + // Check the name for 'PhotolysisRateConstant' + micm::PhotolysisRateConstant* photolysis_rate_const = nullptr; + std::string photolysis_name[] = { "O2_1", "O3_1", "O3_2" }; + + for (short i = 0; i < 3; i++) + { + photolysis_rate_const = dynamic_cast(process_vector[i].rate_constant_.get()); + + EXPECT_TRUE(photolysis_rate_const != nullptr); + EXPECT_EQ(photolysis_rate_const->name_, photolysis_name[i]); + } + + // Check the parameters for 'ArrheniusRateConstant' + micm::ArrheniusRateConstant* arrhenius_rate_const = nullptr; + double A_param[] = { 2.15e-11, 3.3e-11, 8.0e-12, 6.0e-34 }; + double B_param[] = { 0.0, 0.0, 0.0, 2.4 }; + double C_param[] = { 110.0, 55.0, -2060.00, 0.0 }; + double D_param[] = { 300.0, 300.0, 300.0, 300.0 }; + double E_param[] = { 0.0, 0.0, 0.0, 0.0 }; + + for (short i = 3; i < 7; i++) + { + arrhenius_rate_const = dynamic_cast(process_vector[i].rate_constant_.get()); + + EXPECT_TRUE(arrhenius_rate_const != nullptr); + EXPECT_EQ(arrhenius_rate_const->parameters_.A_, A_param[i - 3]); + EXPECT_EQ(arrhenius_rate_const->parameters_.B_, B_param[i - 3]); + EXPECT_EQ(arrhenius_rate_const->parameters_.C_, C_param[i - 3]); + EXPECT_EQ(arrhenius_rate_const->parameters_.D_, D_param[i - 3]); + EXPECT_EQ(arrhenius_rate_const->parameters_.E_, E_param[i - 3]); + } + + // Check the number of custom parameters of 'rate constant' in each 'Process' + std::size_t size_custom_parameters_of_rate_constant_in_each_process[] = { 1, 1, 1, 0, 0, 0, 0 }; + + idx = 0; + std::vector::iterator it; + for (it = solver_params.processes_.begin(); it != solver_params.processes_.end(); it++, idx++) + { + EXPECT_EQ(it->rate_constant_->SizeCustomParameters(), size_custom_parameters_of_rate_constant_in_each_process[idx]); + } +} +#endif \ No newline at end of file diff --git a/test/unit/configure/test_system_builder.cpp b/test/unit/configure/test_system_builder.cpp deleted file mode 100644 index 4dcad6d71..000000000 --- a/test/unit/configure/test_system_builder.cpp +++ /dev/null @@ -1,42 +0,0 @@ -#ifdef USE_JSON -#include -#endif - -#include - -#ifdef USE_JSON -TEST(SystemBuilder, DefaultConstructor){ - micm::SystemBuilder builder{}; -} - -TEST(SystemBuilder, DetectsInvalidConfigFile){ - micm::SystemBuilder builder{}; - EXPECT_ANY_THROW(builder.Build("not_a_config_file.json")); -} - -TEST(SystemBuilder, DetectsInvalidConfigFileAndNoThrowDoesntThrow){ - micm::SystemBuilder builder{}; - EXPECT_NO_THROW(builder.Build("not_a_config_file.json")); - - auto system = builder.Build("not_a_config_file.json"); - ASSERT_TRUE(system.get() == nullptr); -} - -TEST(SystemBuilder, JsonBuilder){ - micm::SystemBuilder builder{}; - auto system = builder.Build("unit_configs/chapman/config.json"); - - EXPECT_TRUE(system != nullptr); - EXPECT_EQ(system->gas_phase_.species_.size(), 9); - for(auto& species : system->gas_phase_.species_){ - if (species.name_ == "M") { - EXPECT_EQ(species.properties_.size(), 0); - } - else{ - EXPECT_EQ(species.properties_.size(), 1); - } - } - EXPECT_EQ(system->photolysis_reactions_.size(), 3); - EXPECT_EQ(system->arrhenius_reactions_.size(), 4); -} -#endif \ No newline at end of file