Skip to content

Commit

Permalink
Surface reaction rate constant (#169)
Browse files Browse the repository at this point in the history
* change species properties to std::map

* add surface reaction rate constant

* add parameters for user defined and surface rate constants

* add cmath include

* add flag for using math constants for windows

* add surface rate config parser

* correcting species assignment operator

* correcting array index

* correcting last array index

---------

Co-authored-by: Kyle Shores <kyle.shores44@gmail.com>
  • Loading branch information
mattldawson and K20shores committed Aug 2, 2023
1 parent 0d9efff commit a690bc1
Show file tree
Hide file tree
Showing 37 changed files with 558 additions and 215 deletions.
109 changes: 72 additions & 37 deletions include/micm/configure/solver_config.hpp
Original file line number Diff line number Diff line change
@@ -1,8 +1,6 @@
/* Copyright (C) 2023 National Center for Atmospheric Research,
*
* SPDX-License-Identifier: Apache-2.0
*
*/
// Copyright (C) 2023 National Center for Atmospheric Research,
//
// SPDX-License-Identifier: Apache-2.0

#pragma once

Expand All @@ -12,13 +10,13 @@
#include <iostream>
#include <micm/process/arrhenius_rate_constant.hpp>
#include <micm/process/branched_rate_constant.hpp>
#include <micm/process/user_defined_rate_constant.hpp>
#include <micm/process/process.hpp>
#include <micm/process/surface_rate_constant.hpp>
#include <micm/process/ternary_chemical_activation_rate_constant.hpp>
#include <micm/process/troe_rate_constant.hpp>
#include <micm/process/tunneling_rate_constant.hpp>
#include <micm/process/user_defined_rate_constant.hpp>
#include <micm/system/phase.hpp>
#include <micm/system/property.hpp>
#include <micm/system/species.hpp>
#include <micm/system/system.hpp>
#include <micm/util/constants.hpp>
Expand Down Expand Up @@ -93,11 +91,10 @@ namespace micm
std::vector<UserDefinedRateConstant> user_defined_rate_arr_;
std::vector<ArrheniusRateConstant> arrhenius_rate_arr_;
std::vector<BranchedRateConstant> branched_rate_arr_;
std::vector<SurfaceRateConstant> surface_rate_arr_;
std::vector<TroeRateConstant> troe_rate_arr_;
std::vector<TernaryChemicalActivationRateConstant> ternary_rate_arr_;
std::vector<TunnelingRateConstant> tunneling_rate_arr_;
std::vector<Species> emission_arr_;
std::vector<Species> first_order_loss_arr_;

// Specific for solver parameters
Phase gas_phase_;
Expand Down Expand Up @@ -282,6 +279,10 @@ namespace micm
{
status = ParseFirstOrderLoss(object);
}
else if (type == "SURFACE")
{
status = ParseSurface(object);
}
else
{
status = ConfigParseStatus::UnknownKey;
Expand All @@ -298,11 +299,6 @@ namespace micm
// required keys
const std::string NAME = "name";

// optional keys
const std::string ABS_TOL = "absolute tolerance";
const std::string MOL_WEIGHT = "molecular weight [kg mol-1]";
const std::string MOL_WEIGHT_UNIT = "kg mol-1";

std::array<std::string, 1> required_keys = { NAME };

// Check if it contains the required key(s)
Expand All @@ -311,24 +307,16 @@ namespace micm
if (!ValidateJsonWithKey(object, key))
return ConfigParseStatus::RequiredKeyNotFound;
}

// Check if it contains optional key(s)
std::string name = object[NAME].get<std::string>();

if (object.contains(ABS_TOL))
{
auto species = Species(name, Property(ABS_TOL, "", object[ABS_TOL].get<double>()));
species_arr_.push_back(species);
}
else if (object.contains(MOL_WEIGHT))
{
auto species = Species(name, Property(MOL_WEIGHT, MOL_WEIGHT_UNIT, object[MOL_WEIGHT].get<double>()));
species_arr_.push_back(species);
}
else
// Load remaining keys as properties
std::map<std::string, double> properties{};
for (auto& [key, value] : object.items())
{
species_arr_.push_back(Species(name));
if (value.is_number_float())
properties[key] = value;
}
species_arr_.push_back(Species(name, properties));

return ConfigParseStatus::Success;
}
Expand Down Expand Up @@ -406,9 +394,10 @@ namespace micm

std::string name = "PHOTO." + object[MUSICA_NAME].get<std::string>();

user_defined_rate_arr_.push_back(UserDefinedRateConstant(name));
user_defined_rate_arr_.push_back(UserDefinedRateConstant({ .label_ = name }));

std::unique_ptr<UserDefinedRateConstant> rate_ptr = std::make_unique<UserDefinedRateConstant>(name);
std::unique_ptr<UserDefinedRateConstant> rate_ptr =
std::make_unique<UserDefinedRateConstant>(UserDefinedRateConstantParameters{ .label_ = name });
processes_.push_back(Process(reactants, products, std::move(rate_ptr), gas_phase_));

return ConfigParseStatus::Success;
Expand Down Expand Up @@ -681,9 +670,10 @@ namespace micm

std::string name = "EMIS." + object[MUSICA_NAME].get<std::string>();

user_defined_rate_arr_.push_back(UserDefinedRateConstant(name));
user_defined_rate_arr_.push_back(UserDefinedRateConstant({ .label_ = name }));

std::unique_ptr<UserDefinedRateConstant> rate_ptr = std::make_unique<UserDefinedRateConstant>(name);
std::unique_ptr<UserDefinedRateConstant> rate_ptr =
std::make_unique<UserDefinedRateConstant>(UserDefinedRateConstantParameters{ .label_ = name });
processes_.push_back(Process(reactants, products, std::move(rate_ptr), gas_phase_));

return ConfigParseStatus::Success;
Expand All @@ -702,17 +692,63 @@ namespace micm
std::string species = object["species"].get<std::string>();
json reactants_object{};
json products_object{};
reactants_object[species] = { { } };
reactants_object[species] = { {} };
auto reactants = ParseReactants(reactants_object);
auto products = ParseProducts(products_object);

std::string name = "LOSS." + object[MUSICA_NAME].get<std::string>();

user_defined_rate_arr_.push_back(UserDefinedRateConstant(name));
user_defined_rate_arr_.push_back(UserDefinedRateConstant({ .label_ = name }));

std::unique_ptr<UserDefinedRateConstant> rate_ptr = std::make_unique<UserDefinedRateConstant>(name);
std::unique_ptr<UserDefinedRateConstant> rate_ptr =
std::make_unique<UserDefinedRateConstant>(UserDefinedRateConstantParameters{ .label_ = name });
processes_.push_back(Process(reactants, products, std::move(rate_ptr), gas_phase_));


return ConfigParseStatus::Success;
}

ConfigParseStatus ParseSurface(const json& object)
{
const std::string REACTANTS = "gas-phase reactant";
const std::string PRODUCTS = "gas-phase products";
const std::string MUSICA_NAME = "MUSICA name";
const std::string PROBABILITY = "reaction probability";
for (const auto& key : { REACTANTS, PRODUCTS, MUSICA_NAME })
{
if (!ValidateJsonWithKey(object, key))
return ConfigParseStatus::RequiredKeyNotFound;
}

std::string species_name = object[REACTANTS].get<std::string>();
json reactants_object{};
reactants_object[species_name] = { {} };
auto reactants = ParseReactants(reactants_object);
auto products = ParseProducts(object[PRODUCTS]);

Species reactant_species = Species("");
for (auto& species : species_arr_)
{
if (species.name_ == species_name)
{
reactant_species = species;
break;
}
}
SurfaceRateConstantParameters parameters{
.label_ = "SURF." + object[MUSICA_NAME].get<std::string>(),
.species_ = reactant_species
};

if (object.contains(PROBABILITY))
{
parameters.reaction_probability_ = object[PROBABILITY].get<double>();
}

surface_rate_arr_.push_back(SurfaceRateConstant(parameters));

std::unique_ptr<SurfaceRateConstant> rate_ptr = std::make_unique<SurfaceRateConstant>(parameters);
processes_.push_back(Process(reactants, products, std::move(rate_ptr), gas_phase_));

return ConfigParseStatus::Success;
}
};
Expand Down Expand Up @@ -748,7 +784,6 @@ namespace micm
return SolverParameters(
std::move(System(std::move(this->gas_phase_), std::move(this->phases_))), std::move(this->processes_));
}

};

} // namespace micm
4 changes: 2 additions & 2 deletions include/micm/process/arrhenius_rate_constant.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@ 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, const std::vector<double>::const_iterator& custom_parameters)
double calculate(const Conditions& conditions, std::vector<double>::const_iterator custom_parameters)
const override;

double calculate(const double& temperature, const double& pressure) const;
Expand All @@ -68,7 +68,7 @@ namespace micm

inline double ArrheniusRateConstant::calculate(
const Conditions& conditions,
const std::vector<double>::const_iterator& custom_parameters) const
std::vector<double>::const_iterator custom_parameters) const
{
return calculate(conditions.temperature_, conditions.pressure_);
}
Expand Down
4 changes: 2 additions & 2 deletions include/micm/process/branched_rate_constant.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -52,7 +52,7 @@ 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, const std::vector<double>::const_iterator& custom_parameters)
double calculate(const Conditions& conditions, std::vector<double>::const_iterator custom_parameters)
const override;

/// @brief Calculate the rate constant
Expand Down Expand Up @@ -88,7 +88,7 @@ namespace micm

inline double BranchedRateConstant::calculate(
const Conditions& conditions,
const std::vector<double>::const_iterator& custom_parameters) const
std::vector<double>::const_iterator custom_parameters) const
{
return calculate(conditions.temperature_, conditions.air_density_);
}
Expand Down
2 changes: 1 addition & 1 deletion include/micm/process/rate_constant.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@ namespace micm
/// @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, const std::vector<double>::const_iterator& custom_parameters)
virtual double calculate(const Conditions& conditions, std::vector<double>::const_iterator custom_parameters)
const
{
return 0;
Expand Down
94 changes: 94 additions & 0 deletions include/micm/process/surface_rate_constant.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,94 @@
// Copyright (C) 2023 National Center for Atmospheric Research,
//
// SPDX-License-Identifier: Apache-2.0
#pragma once

#define _USE_MATH_DEFINES
#include <math.h>

#include <cmath>
#include <micm/process/rate_constant.hpp>
#include <micm/system/species.hpp>
#include <micm/util/constants.hpp>
#include <micm/util/property_keys.hpp>
#include <string>

namespace micm
{
struct SurfaceRateConstantParameters
{
/// @brief Label for the reaction used to identify user-defined parameters
std::string label_;
/// @brief Gas-phase species reacting on surface
Species species_;
/// @brief Reaction probability (0-1) [unitless]
double reaction_probability_{ 1.0 };
};

/// @brief A rate constant for surface reactions
class SurfaceRateConstant : public RateConstant
{
public:
SurfaceRateConstantParameters parameters_;
double diffusion_coefficient_; // [m2 s-1]
double mean_free_speed_factor_; // 8 * gas_constant / ( pi * molecular_weight ) [K-1]

public:
/// @brief
/// @param name A name for this reaction
SurfaceRateConstant(const SurfaceRateConstantParameters& parameters);

/// @brief Deep copy
std::unique_ptr<RateConstant> clone() const override;

/// @brief Returns labels for the surface rate constant parameters:
/// aerosol effective radius [m]
/// aerosol number concentration [# m-3]
/// @return Rate constant labels
std::vector<std::string> CustomParameters() const override;

/// @brief Returns the number of custom parameters
/// @return Number of custom parameters
std::size_t SizeCustomParameters() 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<double>::const_iterator custom_parameters) const override;
};

inline SurfaceRateConstant::SurfaceRateConstant(const SurfaceRateConstantParameters& parameters)
: parameters_(parameters),
diffusion_coefficient_(parameters.species_.properties_.at(GAS_DIFFUSION_COEFFICIENT)),
mean_free_speed_factor_(8.0 * GAS_CONSTANT / (M_PI * parameters.species_.properties_.at(MOLECULAR_WEIGHT)))
{
}

inline std::unique_ptr<RateConstant> SurfaceRateConstant::clone() const
{
return std::unique_ptr<RateConstant>{ new SurfaceRateConstant{ *this } };
}

inline double SurfaceRateConstant::calculate(
const Conditions& conditions,
std::vector<double>::const_iterator custom_parameters) const
{
const double mean_free_speed = std::sqrt(mean_free_speed_factor_ * conditions.temperature_);
const double radius = *(custom_parameters++);
const double number = *(custom_parameters);
return (double)4.0 * number * M_PI * radius * radius /
(radius / diffusion_coefficient_ + 4.0 / (mean_free_speed * parameters_.reaction_probability_));
}

inline std::vector<std::string> SurfaceRateConstant::CustomParameters() const
{
return std::vector<std::string>{ parameters_.label_ + ".effective radius [m]",
parameters_.label_ + ".particle number concentration [# m-3]" };
}

inline std::size_t SurfaceRateConstant::SizeCustomParameters() const
{
return 2;
}
} // namespace micm
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,7 @@ 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, const std::vector<double>::const_iterator& custom_parameters)
double calculate(const Conditions& conditions, std::vector<double>::const_iterator custom_parameters)
const override;

/// @brief Calculate the rate constant
Expand Down Expand Up @@ -78,7 +78,7 @@ namespace micm

inline double TernaryChemicalActivationRateConstant::calculate(
const Conditions& conditions,
const std::vector<double>::const_iterator& custom_parameters) const
std::vector<double>::const_iterator custom_parameters) const
{
return calculate(conditions.temperature_, conditions.air_density_);
}
Expand Down
4 changes: 2 additions & 2 deletions include/micm/process/troe_rate_constant.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,7 @@ 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, const std::vector<double>::const_iterator& custom_parameters)
double calculate(const Conditions& conditions, std::vector<double>::const_iterator custom_parameters)
const override;

/// @brief Calculate the rate constant
Expand All @@ -77,7 +77,7 @@ namespace micm

inline double TroeRateConstant::calculate(
const Conditions& conditions,
const std::vector<double>::const_iterator& custom_parameters) const
std::vector<double>::const_iterator custom_parameters) const
{
return calculate(conditions.temperature_, conditions.air_density_);
}
Expand Down
4 changes: 2 additions & 2 deletions include/micm/process/tunneling_rate_constant.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@ 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, const std::vector<double>::const_iterator& custom_parameters)
double calculate(const Conditions& conditions, std::vector<double>::const_iterator custom_parameters)
const override;

/// @brief Calculate the rate constant
Expand All @@ -66,7 +66,7 @@ namespace micm

inline double TunnelingRateConstant::calculate(
const Conditions& conditions,
const std::vector<double>::const_iterator& custom_parameters) const
std::vector<double>::const_iterator custom_parameters) const
{
return calculate(conditions.temperature_);
}
Expand Down
Loading

0 comments on commit a690bc1

Please sign in to comment.