Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Develop 245 terminator #314

Merged
merged 5 commits into from
Oct 18, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
33 changes: 19 additions & 14 deletions include/micm/process/process.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -36,13 +36,11 @@ namespace micm
/// @param processes The set of processes being solved
/// @param state The solver state to update
template<template<class> class MatrixPolicy>
requires(!VectorizableDense<MatrixPolicy<double>>) static void UpdateState(
const std::vector<Process>& processes,
State<MatrixPolicy>& state);
requires(!VectorizableDense<MatrixPolicy<double>>)
static void UpdateState(const std::vector<Process>& processes, State<MatrixPolicy>& state);
template<template<class> class MatrixPolicy>
requires(VectorizableDense<MatrixPolicy<double>>) static void UpdateState(
const std::vector<Process>& processes,
State<MatrixPolicy>& state);
requires(VectorizableDense<MatrixPolicy<double>>)
static void UpdateState(const std::vector<Process>& processes, State<MatrixPolicy>& state);

friend class ProcessBuilder;
static ProcessBuilder create();
Expand Down Expand Up @@ -92,9 +90,8 @@ namespace micm
};

template<template<class> class MatrixPolicy>
requires(!VectorizableDense<MatrixPolicy<double>>) void Process::UpdateState(
const std::vector<Process>& processes,
State<MatrixPolicy>& state)
requires(!VectorizableDense<MatrixPolicy<double>>)
void Process::UpdateState(const std::vector<Process>& processes, State<MatrixPolicy>& state)
{
for (std::size_t i{}; i < state.custom_rate_parameters_.size(); ++i)
{
Expand All @@ -103,17 +100,20 @@ namespace micm
std::size_t i_rate_constant = 0;
for (auto& process : processes)
{
double fixed_reactants = 1.0;
for (auto& reactant : process.reactants_)
if (reactant.IsParameterized())
fixed_reactants *= reactant.parameterize_(state.conditions_[i]);
state.rate_constants_[i][(i_rate_constant++)] =
process.rate_constant_->calculate(state.conditions_[i], custom_parameters_iter);
process.rate_constant_->calculate(state.conditions_[i], custom_parameters_iter) * fixed_reactants;
custom_parameters_iter += process.rate_constant_->SizeCustomParameters();
}
}
}

template<template<class> class MatrixPolicy>
requires(VectorizableDense<MatrixPolicy<double>>) void Process::UpdateState(
const std::vector<Process>& processes,
State<MatrixPolicy>& state)
requires(VectorizableDense<MatrixPolicy<double>>)
void Process::UpdateState(const std::vector<Process>& processes, State<MatrixPolicy>& state)
{
const auto& v_custom_parameters = state.custom_rate_parameters_.AsVector();
auto& v_rate_constants = state.rate_constants_.AsVector();
Expand All @@ -133,8 +133,13 @@ namespace micm
params[i_param] = v_custom_parameters[offset_params + i_param * L + i_cell];
}
std::vector<double>::const_iterator custom_parameters_iter = params.begin();
double fixed_reactants = 1.0;
for (auto& reactant : process.reactants_)
if (reactant.IsParameterized())
fixed_reactants *= reactant.parameterize_(state.conditions_[i_group * L + i_cell]);
v_rate_constants[offset_rc + i_cell] =
process.rate_constant_->calculate(state.conditions_[i_group * L + i_cell], custom_parameters_iter);
process.rate_constant_->calculate(state.conditions_[i_group * L + i_cell], custom_parameters_iter) *
fixed_reactants;
}
offset_params += params.size() * L;
offset_rc += L;
Expand Down
49 changes: 30 additions & 19 deletions include/micm/process/process_set.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -47,12 +47,13 @@ namespace micm
/// @param state_variables Current state variable values (grid cell, state variable)
/// @param forcing Forcing terms for each state variable (grid cell, state variable)
template<template<class> typename MatrixPolicy>
requires(!VectorizableDense<MatrixPolicy<double>>) void AddForcingTerms(
requires(!VectorizableDense<MatrixPolicy<double>>)
void AddForcingTerms(
const MatrixPolicy<double>& rate_constants,
const MatrixPolicy<double>& state_variables,
MatrixPolicy<double>& forcing) const;
template<template<class> typename MatrixPolicy>
requires VectorizableDense<MatrixPolicy<double>>
requires VectorizableDense<MatrixPolicy<double>>
void AddForcingTerms(
const MatrixPolicy<double>& rate_constants,
const MatrixPolicy<double>& state_variables,
Expand All @@ -63,12 +64,14 @@ namespace micm
/// @param state_variables Current state variable values (grid cell, state variable)
/// @param jacobian Jacobian matrix for the system (grid cell, dependent variable, independent variable)
template<template<class> class MatrixPolicy, template<class> class SparseMatrixPolicy>
requires(!VectorizableDense<MatrixPolicy<double>> || !VectorizableSparse<SparseMatrixPolicy<double>>) void AddJacobianTerms(
requires(!VectorizableDense<MatrixPolicy<double>> || !VectorizableSparse<SparseMatrixPolicy<double>>)
void AddJacobianTerms(
const MatrixPolicy<double>& rate_constants,
const MatrixPolicy<double>& state_variables,
SparseMatrixPolicy<double>& jacobian) const;
template<template<class> class MatrixPolicy, template<class> class SparseMatrixPolicy>
requires(VectorizableDense<MatrixPolicy<double>>&& VectorizableSparse<SparseMatrixPolicy<double>>) void AddJacobianTerms(
requires(VectorizableDense<MatrixPolicy<double>> && VectorizableSparse<SparseMatrixPolicy<double>>)
void AddJacobianTerms(
const MatrixPolicy<double>& rate_constants,
const MatrixPolicy<double>& state_variables,
SparseMatrixPolicy<double>& jacobian) const;
Expand All @@ -84,17 +87,25 @@ namespace micm
{
for (auto& process : processes)
{
number_of_reactants_.push_back(process.reactants_.size());
number_of_products_.push_back(process.products_.size());
std::size_t number_of_reactants = 0;
std::size_t number_of_products = 0;
for (auto& reactant : process.reactants_)
{
if (reactant.IsParameterized())
continue; // Skip reactants that are parameterizations
reactant_ids_.push_back(state.variable_map_.at(reactant.name_));
++number_of_reactants;
}
for (auto& product : process.products_)
{
if (product.first.IsParameterized())
continue; // Skip products that are parameterizations
product_ids_.push_back(state.variable_map_.at(product.first.name_));
yields_.push_back(product.second);
++number_of_products;
}
number_of_reactants_.push_back(number_of_reactants);
number_of_products_.push_back(number_of_products);
}
};

Expand Down Expand Up @@ -147,7 +158,8 @@ namespace micm
}

template<template<class> typename MatrixPolicy>
requires(!VectorizableDense<MatrixPolicy<double>>) inline void ProcessSet::AddForcingTerms(
requires(!VectorizableDense<MatrixPolicy<double>>)
inline void ProcessSet::AddForcingTerms(
const MatrixPolicy<double>& rate_constants,
const MatrixPolicy<double>& state_variables,
MatrixPolicy<double>& forcing) const
Expand Down Expand Up @@ -183,7 +195,7 @@ namespace micm
};

template<template<class> typename MatrixPolicy>
requires VectorizableDense<MatrixPolicy<double>>
requires VectorizableDense<MatrixPolicy<double>>
inline void ProcessSet::AddForcingTerms(
const MatrixPolicy<double>& rate_constants,
const MatrixPolicy<double>& state_variables,
Expand Down Expand Up @@ -224,12 +236,11 @@ namespace micm
}

template<template<class> class MatrixPolicy, template<class> class SparseMatrixPolicy>
requires(
!VectorizableDense<MatrixPolicy<double>> || !VectorizableSparse<SparseMatrixPolicy<double>>) inline void ProcessSet::
AddJacobianTerms(
const MatrixPolicy<double>& rate_constants,
const MatrixPolicy<double>& state_variables,
SparseMatrixPolicy<double>& jacobian) const
requires(!VectorizableDense<MatrixPolicy<double>> || !VectorizableSparse<SparseMatrixPolicy<double>>)
inline void ProcessSet::AddJacobianTerms(
const MatrixPolicy<double>& rate_constants,
const MatrixPolicy<double>& state_variables,
SparseMatrixPolicy<double>& jacobian) const
{
auto cell_jacobian = jacobian.AsVector().begin();

Expand Down Expand Up @@ -271,11 +282,11 @@ namespace micm
}

template<template<class> class MatrixPolicy, template<class> class SparseMatrixPolicy>
requires(VectorizableDense<MatrixPolicy<double>>&& VectorizableSparse<SparseMatrixPolicy<double>>) inline void ProcessSet::
AddJacobianTerms(
const MatrixPolicy<double>& rate_constants,
const MatrixPolicy<double>& state_variables,
SparseMatrixPolicy<double>& jacobian) const
requires(VectorizableDense<MatrixPolicy<double>> && VectorizableSparse<SparseMatrixPolicy<double>>)
inline void ProcessSet::AddJacobianTerms(
const MatrixPolicy<double>& rate_constants,
const MatrixPolicy<double>& state_variables,
SparseMatrixPolicy<double>& jacobian) const
{
const auto& v_rate_constants = rate_constants.AsVector();
const auto& v_state_variables = state_variables.AsVector();
Expand Down
11 changes: 7 additions & 4 deletions include/micm/solver/rosenbrock.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -50,9 +50,11 @@ namespace micm
double rejection_factor_decrease_{ 0.1 }; // used to decrease the step after 2 successive rejections
double safety_factor_{ 0.9 }; // safety factor in new step size computation

double h_min_{ 0 }; // step size min
double h_max_{ 0.5 }; // step size max
double h_start_{ 0.005 }; // step size start
double h_min_{ 0.0 }; // step size min [s]
double h_max_{
0.0
}; // step size max [s] (if zero or greater than the solver time-step, the time-step passed to the solver will be used)
double h_start_{ 0.0 }; // step size start [s] (if zero, 1.0e-6 will be used, if greater than h_max, h_max will be used)

// Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
// or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
Expand Down Expand Up @@ -173,7 +175,8 @@ namespace micm
uint64_t decompositions{};
/// @brief The number of linear solvers
uint64_t solves{};
/// @brief The number of times a singular matrix is detected. For now, this will always be zero as we assume the matrix is never singular
/// @brief The number of times a singular matrix is detected. For now, this will always be zero as we assume the matrix
/// is never singular
uint64_t singular{};
/// @brief The cumulative amount of time spent calculating the forcing function
std::chrono::duration<double, std::nano> total_forcing_time{};
Expand Down
39 changes: 22 additions & 17 deletions include/micm/solver/rosenbrock.inl
Original file line number Diff line number Diff line change
@@ -1,17 +1,20 @@
// Copyright (C) 2023 National Center for Atmospheric Research
// SPDX-License-Identifier: Apache-2.0

#define TIMED_METHOD(assigned_increment, time_it, method, ...) \
{ \
if constexpr (time_it) { \
auto start = std::chrono::high_resolution_clock::now(); \
method(__VA_ARGS__); \
auto end = std::chrono::high_resolution_clock::now(); \
assigned_increment += std::chrono::duration_cast<std::chrono::nanoseconds>(end - start); \
} else { \
method(__VA_ARGS__); \
} \
} \
#define TIMED_METHOD(assigned_increment, time_it, method, ...) \
{ \
if constexpr (time_it) \
{ \
auto start = std::chrono::high_resolution_clock::now(); \
method(__VA_ARGS__); \
auto end = std::chrono::high_resolution_clock::now(); \
assigned_increment += std::chrono::duration_cast<std::chrono::nanoseconds>(end - start); \
} \
else \
{ \
method(__VA_ARGS__); \
} \
}

namespace micm
{
Expand Down Expand Up @@ -517,9 +520,9 @@ namespace micm
MatrixPolicy<double> forcing(Y.size(), Y[0].size(), 0.0);
MatrixPolicy<double> temp(Y.size(), Y[0].size(), 0.0);
std::vector<MatrixPolicy<double>> K{};

parameters_.h_max_ = time_step;
parameters_.h_start_ = std::max(parameters_.h_min_, delta_min_);
const double h_max = parameters_.h_max_ == 0.0 ? time_step : std::min(time_step, parameters_.h_max_);
const double h_start =
parameters_.h_start_ == 0.0 ? std::max(parameters_.h_min_, delta_min_) : std::min(h_max, parameters_.h_start_);

stats_.Reset();
UpdateState(state);
Expand All @@ -528,12 +531,14 @@ namespace micm
K.push_back(MatrixPolicy<double>(Y.size(), Y[0].size(), 0.0));

double present_time = 0.0;
double H =
std::min(std::max(std::abs(parameters_.h_min_), std::abs(parameters_.h_start_)), std::abs(parameters_.h_max_));
double H = std::min(std::max(std::abs(parameters_.h_min_), std::abs(h_start)), std::abs(h_max));

if (std::abs(H) <= 10 * parameters_.round_off_)
H = delta_min_;

// TODO: the logic above this point should be moved to the constructor and should return an error
// if the parameters are invalid (e.g., h_min > h_max or h_start > h_max)

bool reject_last_h = false;
bool reject_more_h = false;

Expand Down Expand Up @@ -639,7 +644,7 @@ namespace micm
stats_.accepted += 1;
present_time = present_time + H;
Y.AsVector().assign(Ynew.AsVector().begin(), Ynew.AsVector().end());
Hnew = std::max(parameters_.h_min_, std::min(Hnew, parameters_.h_max_));
Hnew = std::max(parameters_.h_min_, std::min(Hnew, h_max));
if (reject_last_h)
{
// No step size increase after a rejected step
Expand Down
16 changes: 16 additions & 0 deletions include/micm/system/phase.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -49,5 +49,21 @@ namespace micm
species_ = other.species_;
return *this;
}

/// @brief Returns the number of non-parameterized species
std::size_t StateSize() const
{
return std::count_if(species_.begin(), species_.end(), [](const Species& s) { return !s.IsParameterized(); });
}

/// @brief Returns a set of unique names for each non-parameterized species
std::vector<std::string> UniqueNames() const
{
std::vector<std::string> names{};
for (const auto& species : species_)
if (!species.IsParameterized())
names.push_back(species.name_);
return names;
}
};
} // namespace micm
30 changes: 29 additions & 1 deletion include/micm/system/species.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -3,9 +3,12 @@
// SPDX-License-Identifier: Apache-2.0
#pragma once

#include <functional>
#include <string>
#include <vector>

#include "conditions.hpp"

namespace micm
{

Expand All @@ -19,6 +22,12 @@ namespace micm
/// @brief A list of properties of this species
std::map<std::string, double> properties_;

/// @brief A function that if provided will be used to parameterize
/// the concentration of this species during solving.
/// Species with this function defined will be excluded from
/// the solver state.
std::function<double(const Conditions)> parameterize_{ nullptr };

/// @brief Copy assignment
/// @param other species to copy
Species& operator=(const Species& other);
Expand All @@ -35,6 +44,12 @@ namespace micm
/// @param name The name of the species
/// @param properties The properties of the species
Species(const std::string& name, const std::map<std::string, double>& properties);

/// @brief Returns whether a species is parameterized
bool IsParameterized() const;

/// @brief Return a Species instance parameterized on air density
static Species ThirdBody();
};

inline Species& Species::operator=(const Species& other)
Expand All @@ -46,13 +61,15 @@ namespace micm

name_ = other.name_;
properties_ = other.properties_; // This performs a shallow copy
parameterize_ = other.parameterize_;

return *this;
}

inline Species::Species(const Species& other)
: name_(other.name_),
properties_(other.properties_){};
properties_(other.properties_),
parameterize_(other.parameterize_){};

inline Species::Species(const std::string& name)
: name_(name){};
Expand All @@ -61,4 +78,15 @@ namespace micm
: name_(name),
properties_(properties){};

inline bool Species::IsParameterized() const
{
return parameterize_ != nullptr;
}

inline Species Species::ThirdBody()
{
Species third_body{ "M" };
third_body.parameterize_ = [](const Conditions& c) { return c.air_density_; };
return third_body;
}
} // namespace micm
Loading