diff --git a/include/micm/solver/linear_solver.hpp b/include/micm/solver/linear_solver.hpp index 492998c08..61a399e68 100644 --- a/include/micm/solver/linear_solver.hpp +++ b/include/micm/solver/linear_solver.hpp @@ -70,220 +70,6 @@ namespace micm MatrixPolicy& x); }; - template class MatrixPolicy> - std::vector DiagonalMarkowitzReorder(const MatrixPolicy& matrix) - { - const std::size_t order = matrix.size(); - std::vector perm(order); - for (std::size_t i = 0; i < order; ++i) - perm[i] = i; - MatrixPolicy pattern = matrix; - for (std::size_t row = 0; row < (order - 1); ++row) - { - std::size_t beta = std::pow((order - 1), 2); - std::size_t max_row = row; - for (std::size_t col = row; col < order; ++col) - { - std::size_t count_a = 0; - std::size_t count_b = 0; - for (std::size_t i = row; i < order; ++i) - { - count_a += (pattern[col][i] == 0 ? 0 : 1); - count_b += (pattern[i][col] == 0 ? 0 : 1); - } - std::size_t count = (count_a - 1) * (count_b - 1); - if (count < beta) - { - beta = count; - max_row = col; - } - } - // Swap row and max_row - if (max_row != row) - { - for (std::size_t i = row; i < order; ++i) - std::swap(pattern[row][i], pattern[max_row][i]); - for (std::size_t i = row; i < order; ++i) - std::swap(pattern[i][row], pattern[i][max_row]); - std::swap(perm[row], perm[max_row]); - } - for (std::size_t col = row + 1; col < order; ++col) - if (pattern[row][col]) - for (std::size_t i = row + 1; i < order; ++i) - pattern[i][col] = pattern[i][row] || pattern[i][col]; - } - return perm; - } - - template class SparseMatrixPolicy> - inline LinearSolver::LinearSolver(const SparseMatrixPolicy& matrix, T initial_value) - : nLij_Lii_(), - Lij_yj_(), - nUij_Uii_(), - Uij_xj_(), - lu_decomp_(matrix) - { - auto lu = lu_decomp_.GetLUMatrices(matrix, initial_value); - lower_matrix_ = std::move(lu.first); - upper_matrix_ = std::move(lu.second); - for (std::size_t i = 0; i < lower_matrix_[0].size(); ++i) - { - std::size_t nLij = 0; - for (std::size_t j_id = lower_matrix_.RowStartVector()[i]; j_id < lower_matrix_.RowStartVector()[i + 1]; ++j_id) - { - std::size_t j = lower_matrix_.RowIdsVector()[j_id]; - if (j >= i) - break; - Lij_yj_.push_back(std::make_pair(lower_matrix_.VectorIndex(0, i, j), j)); - ++nLij; - } - // There must always be a non-zero element on the diagonal - nLij_Lii_.push_back(std::make_pair(nLij, lower_matrix_.VectorIndex(0, i, i))); - } - for (std::size_t i = upper_matrix_[0].size() - 1; i != static_cast(-1); --i) - { - std::size_t nUij = 0; - for (std::size_t j_id = upper_matrix_.RowStartVector()[i]; j_id < upper_matrix_.RowStartVector()[i + 1]; ++j_id) - { - std::size_t j = upper_matrix_.RowIdsVector()[j_id]; - if (j <= i) - continue; - Uij_xj_.push_back(std::make_pair(upper_matrix_.VectorIndex(0, i, j), j)); - ++nUij; - } - // There must always be a non-zero element on the diagonal - nUij_Uii_.push_back(std::make_pair(nUij, upper_matrix_.VectorIndex(0, i, i))); - } - }; - - template class SparseMatrixPolicy> - inline void LinearSolver::Factor(SparseMatrixPolicy& matrix) - { - lu_decomp_.Decompose(matrix, lower_matrix_, upper_matrix_); - } - - template class SparseMatrixPolicy> - template class MatrixPolicy> - requires(!VectorizableDense> || !VectorizableSparse>) inline void LinearSolver< - T, - SparseMatrixPolicy>::Solve(const MatrixPolicy& b, MatrixPolicy& x) - { - for (std::size_t i_cell = 0; i_cell < b.size(); ++i_cell) - { - auto b_cell = b[i_cell]; - auto x_cell = x[i_cell]; - const std::size_t lower_grid_offset = i_cell * lower_matrix_.FlatBlockSize(); - const std::size_t upper_grid_offset = i_cell * upper_matrix_.FlatBlockSize(); - auto& y_cell = x_cell; // Alias x for consistency with equations, but to reuse memory - { - auto b_elem = b_cell.begin(); - auto y_elem = y_cell.begin(); - auto Lij_yj = Lij_yj_.begin(); - for (auto& nLij_Lii : nLij_Lii_) - { - *y_elem = *(b_elem++); - for (std::size_t i = 0; i < nLij_Lii.first; ++i) - { - *y_elem -= lower_matrix_.AsVector()[lower_grid_offset + (*Lij_yj).first] * y_cell[(*Lij_yj).second]; - ++Lij_yj; - } - *(y_elem++) /= lower_matrix_.AsVector()[lower_grid_offset + nLij_Lii.second]; - } - } - { - auto y_elem = std::next(y_cell.end(), -1); - auto x_elem = std::next(x_cell.end(), -1); - auto Uij_xj = Uij_xj_.begin(); - for (auto& nUij_Uii : nUij_Uii_) - { - *x_elem = *(y_elem); - // don't iterate before the beginning of the vector - if (y_elem != y_cell.begin()) - { - --y_elem; - } - - for (std::size_t i = 0; i < nUij_Uii.first; ++i) - { - *x_elem -= upper_matrix_.AsVector()[upper_grid_offset + (*Uij_xj).first] * x_cell[(*Uij_xj).second]; - ++Uij_xj; - } - - // don't iterate before the beginning of the vector - *(x_elem) /= upper_matrix_.AsVector()[upper_grid_offset + nUij_Uii.second]; - if (x_elem != x_cell.begin()) - { - --x_elem; - } - } - } - } - } - - template class SparseMatrixPolicy> - template class MatrixPolicy> - requires(VectorizableDense>&& VectorizableSparse>) inline void LinearSolver::Solve(const MatrixPolicy& b, MatrixPolicy& x) - { - const std::size_t n_cells = b.GroupVectorSize(); - // Loop over groups of blocks - for (std::size_t i_group = 0; i_group < b.NumberOfGroups(); ++i_group) - { - auto b_group = std::next(b.AsVector().begin(), i_group * b.GroupSize()); - auto x_group = std::next(x.AsVector().begin(), i_group * x.GroupSize()); - auto L_group = - std::next(lower_matrix_.AsVector().begin(), i_group * lower_matrix_.GroupSize(lower_matrix_.FlatBlockSize())); - auto U_group = - std::next(upper_matrix_.AsVector().begin(), i_group * upper_matrix_.GroupSize(upper_matrix_.FlatBlockSize())); - auto y_group = x_group; // Alias x for consistency with equations, but to reuse memory - { - auto b_elem = b_group; - auto y_elem = y_group; - auto Lij_yj = Lij_yj_.begin(); - for (auto& nLij_Lii : nLij_Lii_) - { - for (std::size_t i_cell = 0; i_cell < n_cells; ++i_cell) - y_elem[i_cell] = b_elem[i_cell]; - b_elem += n_cells; - for (std::size_t i = 0; i < nLij_Lii.first; ++i) - { - for (std::size_t i_cell = 0; i_cell < n_cells; ++i_cell) - y_elem[i_cell] -= L_group[(*Lij_yj).first + i_cell] * y_group[(*Lij_yj).second * n_cells + i_cell]; - ++Lij_yj; - } - for (std::size_t i_cell = 0; i_cell < n_cells; ++i_cell) - y_elem[i_cell] /= L_group[nLij_Lii.second + i_cell]; - y_elem += n_cells; - } - } - { - auto y_elem = std::next(y_group, x.GroupSize() - n_cells); - auto x_elem = std::next(x_group, x.GroupSize() - n_cells); - auto Uij_xj = Uij_xj_.begin(); - for (auto& nUij_Uii : nUij_Uii_) - { - for (std::size_t i_cell = 0; i_cell < n_cells; ++i_cell) - x_elem[i_cell] = y_elem[i_cell]; - - // don't iterate before the beginning of the vector - std::size_t y_elem_distance = std::distance(x.AsVector().begin(), y_elem); - y_elem -= std::min(n_cells, y_elem_distance); - - for (std::size_t i = 0; i < nUij_Uii.first; ++i) - { - for (std::size_t i_cell = 0; i_cell < n_cells; ++i_cell) - x_elem[i_cell] -= U_group[(*Uij_xj).first + i_cell] * x_group[(*Uij_xj).second * n_cells + i_cell]; - ++Uij_xj; - } - for (std::size_t i_cell = 0; i_cell < n_cells; ++i_cell) - x_elem[i_cell] /= U_group[nUij_Uii.second + i_cell]; - - // don't iterate before the beginning of the vector - std::size_t x_elem_distance = std::distance(x.AsVector().begin(), x_elem); - x_elem -= std::min(n_cells, x_elem_distance); - } - } - } - } +} // namespace micm -} // namespace micm \ No newline at end of file +#include "linear_solver.inl" \ No newline at end of file diff --git a/include/micm/solver/linear_solver.inl b/include/micm/solver/linear_solver.inl new file mode 100644 index 000000000..7ca9ff7c7 --- /dev/null +++ b/include/micm/solver/linear_solver.inl @@ -0,0 +1,220 @@ +namespace micm +{ + + template class MatrixPolicy> + inline std::vector DiagonalMarkowitzReorder(const MatrixPolicy& matrix) + { + const std::size_t order = matrix.size(); + std::vector perm(order); + for (std::size_t i = 0; i < order; ++i) + perm[i] = i; + MatrixPolicy pattern = matrix; + for (std::size_t row = 0; row < (order - 1); ++row) + { + std::size_t beta = std::pow((order - 1), 2); + std::size_t max_row = row; + for (std::size_t col = row; col < order; ++col) + { + std::size_t count_a = 0; + std::size_t count_b = 0; + for (std::size_t i = row; i < order; ++i) + { + count_a += (pattern[col][i] == 0 ? 0 : 1); + count_b += (pattern[i][col] == 0 ? 0 : 1); + } + std::size_t count = (count_a - 1) * (count_b - 1); + if (count < beta) + { + beta = count; + max_row = col; + } + } + // Swap row and max_row + if (max_row != row) + { + for (std::size_t i = row; i < order; ++i) + std::swap(pattern[row][i], pattern[max_row][i]); + for (std::size_t i = row; i < order; ++i) + std::swap(pattern[i][row], pattern[i][max_row]); + std::swap(perm[row], perm[max_row]); + } + for (std::size_t col = row + 1; col < order; ++col) + if (pattern[row][col]) + for (std::size_t i = row + 1; i < order; ++i) + pattern[i][col] = pattern[i][row] || pattern[i][col]; + } + return perm; + } + + template class SparseMatrixPolicy> + inline LinearSolver::LinearSolver(const SparseMatrixPolicy& matrix, T initial_value) + : nLij_Lii_(), + Lij_yj_(), + nUij_Uii_(), + Uij_xj_(), + lu_decomp_(matrix) + { + auto lu = lu_decomp_.GetLUMatrices(matrix, initial_value); + lower_matrix_ = std::move(lu.first); + upper_matrix_ = std::move(lu.second); + for (std::size_t i = 0; i < lower_matrix_[0].size(); ++i) + { + std::size_t nLij = 0; + for (std::size_t j_id = lower_matrix_.RowStartVector()[i]; j_id < lower_matrix_.RowStartVector()[i + 1]; ++j_id) + { + std::size_t j = lower_matrix_.RowIdsVector()[j_id]; + if (j >= i) + break; + Lij_yj_.push_back(std::make_pair(lower_matrix_.VectorIndex(0, i, j), j)); + ++nLij; + } + // There must always be a non-zero element on the diagonal + nLij_Lii_.push_back(std::make_pair(nLij, lower_matrix_.VectorIndex(0, i, i))); + } + for (std::size_t i = upper_matrix_[0].size() - 1; i != static_cast(-1); --i) + { + std::size_t nUij = 0; + for (std::size_t j_id = upper_matrix_.RowStartVector()[i]; j_id < upper_matrix_.RowStartVector()[i + 1]; ++j_id) + { + std::size_t j = upper_matrix_.RowIdsVector()[j_id]; + if (j <= i) + continue; + Uij_xj_.push_back(std::make_pair(upper_matrix_.VectorIndex(0, i, j), j)); + ++nUij; + } + // There must always be a non-zero element on the diagonal + nUij_Uii_.push_back(std::make_pair(nUij, upper_matrix_.VectorIndex(0, i, i))); + } + }; + + template class SparseMatrixPolicy> + inline void LinearSolver::Factor(SparseMatrixPolicy& matrix) + { + lu_decomp_.Decompose(matrix, lower_matrix_, upper_matrix_); + } + + template class SparseMatrixPolicy> + template class MatrixPolicy> + requires(!VectorizableDense> || !VectorizableSparse>) inline void LinearSolver< + T, + SparseMatrixPolicy>::Solve(const MatrixPolicy& b, MatrixPolicy& x) + { + for (std::size_t i_cell = 0; i_cell < b.size(); ++i_cell) + { + auto b_cell = b[i_cell]; + auto x_cell = x[i_cell]; + const std::size_t lower_grid_offset = i_cell * lower_matrix_.FlatBlockSize(); + const std::size_t upper_grid_offset = i_cell * upper_matrix_.FlatBlockSize(); + auto& y_cell = x_cell; // Alias x for consistency with equations, but to reuse memory + { + auto b_elem = b_cell.begin(); + auto y_elem = y_cell.begin(); + auto Lij_yj = Lij_yj_.begin(); + for (auto& nLij_Lii : nLij_Lii_) + { + *y_elem = *(b_elem++); + for (std::size_t i = 0; i < nLij_Lii.first; ++i) + { + *y_elem -= lower_matrix_.AsVector()[lower_grid_offset + (*Lij_yj).first] * y_cell[(*Lij_yj).second]; + ++Lij_yj; + } + *(y_elem++) /= lower_matrix_.AsVector()[lower_grid_offset + nLij_Lii.second]; + } + } + { + auto y_elem = std::next(y_cell.end(), -1); + auto x_elem = std::next(x_cell.end(), -1); + auto Uij_xj = Uij_xj_.begin(); + for (auto& nUij_Uii : nUij_Uii_) + { + *x_elem = *(y_elem); + // don't iterate before the beginning of the vector + if (y_elem != y_cell.begin()) + { + --y_elem; + } + + for (std::size_t i = 0; i < nUij_Uii.first; ++i) + { + *x_elem -= upper_matrix_.AsVector()[upper_grid_offset + (*Uij_xj).first] * x_cell[(*Uij_xj).second]; + ++Uij_xj; + } + + // don't iterate before the beginning of the vector + *(x_elem) /= upper_matrix_.AsVector()[upper_grid_offset + nUij_Uii.second]; + if (x_elem != x_cell.begin()) + { + --x_elem; + } + } + } + } + } + + template class SparseMatrixPolicy> + template class MatrixPolicy> + requires(VectorizableDense>&& VectorizableSparse>) inline void LinearSolver::Solve(const MatrixPolicy& b, MatrixPolicy& x) + { + const std::size_t n_cells = b.GroupVectorSize(); + // Loop over groups of blocks + for (std::size_t i_group = 0; i_group < b.NumberOfGroups(); ++i_group) + { + auto b_group = std::next(b.AsVector().begin(), i_group * b.GroupSize()); + auto x_group = std::next(x.AsVector().begin(), i_group * x.GroupSize()); + auto L_group = + std::next(lower_matrix_.AsVector().begin(), i_group * lower_matrix_.GroupSize(lower_matrix_.FlatBlockSize())); + auto U_group = + std::next(upper_matrix_.AsVector().begin(), i_group * upper_matrix_.GroupSize(upper_matrix_.FlatBlockSize())); + auto y_group = x_group; // Alias x for consistency with equations, but to reuse memory + { + auto b_elem = b_group; + auto y_elem = y_group; + auto Lij_yj = Lij_yj_.begin(); + for (auto& nLij_Lii : nLij_Lii_) + { + for (std::size_t i_cell = 0; i_cell < n_cells; ++i_cell) + y_elem[i_cell] = b_elem[i_cell]; + b_elem += n_cells; + for (std::size_t i = 0; i < nLij_Lii.first; ++i) + { + for (std::size_t i_cell = 0; i_cell < n_cells; ++i_cell) + y_elem[i_cell] -= L_group[(*Lij_yj).first + i_cell] * y_group[(*Lij_yj).second * n_cells + i_cell]; + ++Lij_yj; + } + for (std::size_t i_cell = 0; i_cell < n_cells; ++i_cell) + y_elem[i_cell] /= L_group[nLij_Lii.second + i_cell]; + y_elem += n_cells; + } + } + { + auto y_elem = std::next(y_group, x.GroupSize() - n_cells); + auto x_elem = std::next(x_group, x.GroupSize() - n_cells); + auto Uij_xj = Uij_xj_.begin(); + for (auto& nUij_Uii : nUij_Uii_) + { + for (std::size_t i_cell = 0; i_cell < n_cells; ++i_cell) + x_elem[i_cell] = y_elem[i_cell]; + + // don't iterate before the beginning of the vector + std::size_t y_elem_distance = std::distance(x.AsVector().begin(), y_elem); + y_elem -= std::min(n_cells, y_elem_distance); + + for (std::size_t i = 0; i < nUij_Uii.first; ++i) + { + for (std::size_t i_cell = 0; i_cell < n_cells; ++i_cell) + x_elem[i_cell] -= U_group[(*Uij_xj).first + i_cell] * x_group[(*Uij_xj).second * n_cells + i_cell]; + ++Uij_xj; + } + for (std::size_t i_cell = 0; i_cell < n_cells; ++i_cell) + x_elem[i_cell] /= U_group[nUij_Uii.second + i_cell]; + + // don't iterate before the beginning of the vector + std::size_t x_elem_distance = std::distance(x.AsVector().begin(), x_elem); + x_elem -= std::min(n_cells, x_elem_distance); + } + } + } + } + +} // namespace micm \ No newline at end of file diff --git a/include/micm/solver/lu_decomposition.hpp b/include/micm/solver/lu_decomposition.hpp index 8866994ba..dea4cbfde 100644 --- a/include/micm/solver/lu_decomposition.hpp +++ b/include/micm/solver/lu_decomposition.hpp @@ -101,265 +101,6 @@ namespace micm SparseMatrixPolicy& U) const; }; - inline LuDecomposition::LuDecomposition() - { - } - - template - inline LuDecomposition::LuDecomposition(const SparseMatrix& matrix) - { - std::size_t n = matrix[0].size(); - auto LU = GetLUMatrices(matrix, T{}); - const auto& L_row_start = LU.first.RowStartVector(); - const auto& L_row_ids = LU.first.RowIdsVector(); - const auto& U_row_start = LU.second.RowStartVector(); - const auto& U_row_ids = LU.second.RowIdsVector(); - for (std::size_t i = 0; i < matrix[0].size(); ++i) - { - std::pair iLU(0, 0); - // Upper triangular matrix - for (std::size_t k = i; k < n; ++k) - { - std::size_t nkj = 0; - for (std::size_t j_id = L_row_start[i]; j_id < L_row_start[i + 1]; ++j_id) - { - std::size_t j = L_row_ids[j_id]; - if (j >= i) - break; - if (LU.second.IsZero(j, k)) - continue; - ++nkj; - lij_ujk_.push_back(std::make_pair(LU.first.VectorIndex(0, i, j), LU.second.VectorIndex(0, j, k))); - } - if (matrix.IsZero(i, k)) - { - if (nkj == 0 && k != i) - continue; - do_aik_.push_back(false); - } - else - { - do_aik_.push_back(true); - aik_.push_back(matrix.VectorIndex(0, i, k)); - } - uik_nkj_.push_back(std::make_pair(LU.second.VectorIndex(0, i, k), nkj)); - ++(iLU.second); - } - // Lower triangular matrix - lki_nkj_.push_back(std::make_pair(LU.first.VectorIndex(0, i, i), 0)); - for (std::size_t k = i + 1; k < n; ++k) - { - std::size_t nkj = 0; - for (std::size_t j_id = L_row_start[k]; j_id < L_row_start[k + 1]; ++j_id) - { - std::size_t j = L_row_ids[j_id]; - if (j >= i) - break; - if (LU.second.IsZero(j, i)) - continue; - ++nkj; - lkj_uji_.push_back(std::make_pair(LU.first.VectorIndex(0, k, j), LU.second.VectorIndex(0, j, i))); - } - if (matrix.IsZero(k, i)) - { - if (nkj == 0) - continue; - do_aki_.push_back(false); - } - else - { - do_aki_.push_back(true); - aki_.push_back(matrix.VectorIndex(0, k, i)); - } - uii_.push_back(LU.second.VectorIndex(0, i, i)); - lki_nkj_.push_back(std::make_pair(LU.first.VectorIndex(0, k, i), nkj)); - ++(iLU.first); - } - niLU_.push_back(iLU); - } - } - - template - inline std::pair, SparseMatrix> LuDecomposition::GetLUMatrices( - const SparseMatrix& A, - T initial_value) - { - std::size_t n = A[0].size(); - std::set> L_ids, U_ids; - const auto& row_start = A.RowStartVector(); - const auto& row_ids = A.RowIdsVector(); - for (std::size_t i = 0; i < n; ++i) - { - // Upper triangular matrix - for (std::size_t k = i; k < n; ++k) - { - if (!A.IsZero(i, k) || k == i) - { - U_ids.insert(std::make_pair(i, k)); - continue; - } - for (std::size_t j = 0; j < i; ++j) - { - if (L_ids.find(std::make_pair(i, j)) != L_ids.end() && U_ids.find(std::make_pair(j, k)) != U_ids.end()) - { - U_ids.insert(std::make_pair(i, k)); - break; - } - } - } - // Lower triangular matrix - for (std::size_t k = i; k < n; ++k) - { - if (!A.IsZero(k, i) || k == i) - { - L_ids.insert(std::make_pair(k, i)); - continue; - } - for (std::size_t j = 0; j < i; ++j) - { - if (L_ids.find(std::make_pair(k, j)) != L_ids.end() && U_ids.find(std::make_pair(j, i)) != U_ids.end()) - { - L_ids.insert(std::make_pair(k, i)); - break; - } - } - } - } - auto L_builder = - micm::SparseMatrix::create(n).number_of_blocks(A.size()).initial_value(initial_value); - for (auto& pair : L_ids) - { - L_builder = L_builder.with_element(pair.first, pair.second); - } - auto U_builder = - micm::SparseMatrix::create(n).number_of_blocks(A.size()).initial_value(initial_value); - for (auto& pair : U_ids) - { - U_builder = U_builder.with_element(pair.first, pair.second); - } - std::pair, SparseMatrix> LU(L_builder, U_builder); - return LU; - } - - template class SparseMatrixPolicy> - requires(!VectorizableSparse>) void LuDecomposition::Decompose( - const SparseMatrixPolicy& A, - SparseMatrixPolicy& L, - SparseMatrixPolicy& U) const - { - // Loop over blocks - for (std::size_t i_block = 0; i_block < A.size(); ++i_block) - { - auto A_vector = std::next(A.AsVector().begin(), i_block * A.FlatBlockSize()); - auto L_vector = std::next(L.AsVector().begin(), i_block * L.FlatBlockSize()); - auto U_vector = std::next(U.AsVector().begin(), i_block * U.FlatBlockSize()); - auto do_aik = do_aik_.begin(); - auto aik = aik_.begin(); - auto uik_nkj = uik_nkj_.begin(); - auto lij_ujk = lij_ujk_.begin(); - auto do_aki = do_aki_.begin(); - auto aki = aki_.begin(); - auto lki_nkj = lki_nkj_.begin(); - auto lkj_uji = lkj_uji_.begin(); - auto uii = uii_.begin(); - for (auto& inLU : niLU_) - { - // Upper trianglur matrix - for (std::size_t iU = 0; iU < inLU.second; ++iU) - { - if (*(do_aik++)) - U_vector[uik_nkj->first] = A_vector[*(aik++)]; - for (std::size_t ikj = 0; ikj < uik_nkj->second; ++ikj) - { - U_vector[uik_nkj->first] -= L_vector[lij_ujk->first] * U_vector[lij_ujk->second]; - ++lij_ujk; - } - ++uik_nkj; - } - // Lower triangular matrix - L_vector[(lki_nkj++)->first] = 1.0; - for (std::size_t iL = 0; iL < inLU.first; ++iL) - { - if (*(do_aki++)) - L_vector[lki_nkj->first] = A_vector[*(aki++)]; - for (std::size_t ikj = 0; ikj < lki_nkj->second; ++ikj) - { - L_vector[lki_nkj->first] -= L_vector[lkj_uji->first] * U_vector[lkj_uji->second]; - ++lkj_uji; - } - L_vector[lki_nkj->first] /= U_vector[*uii]; - ++lki_nkj; - ++uii; - } - } - } - } - - template class SparseMatrixPolicy> - requires(VectorizableSparse>) void LuDecomposition::Decompose( - const SparseMatrixPolicy& A, - SparseMatrixPolicy& L, - SparseMatrixPolicy& U) const - { - const std::size_t n_cells = A.GroupVectorSize(); - // Loop over groups of blocks - for (std::size_t i_group = 0; i_group < A.NumberOfGroups(A.size()); ++i_group) - { - auto A_vector = std::next(A.AsVector().begin(), i_group * A.GroupSize(A.FlatBlockSize())); - auto L_vector = std::next(L.AsVector().begin(), i_group * L.GroupSize(L.FlatBlockSize())); - auto U_vector = std::next(U.AsVector().begin(), i_group * U.GroupSize(U.FlatBlockSize())); - auto do_aik = do_aik_.begin(); - auto aik = aik_.begin(); - auto uik_nkj = uik_nkj_.begin(); - auto lij_ujk = lij_ujk_.begin(); - auto do_aki = do_aki_.begin(); - auto aki = aki_.begin(); - auto lki_nkj = lki_nkj_.begin(); - auto lkj_uji = lkj_uji_.begin(); - auto uii = uii_.begin(); - for (auto& inLU : niLU_) - { - // Upper trianglur matrix - for (std::size_t iU = 0; iU < inLU.second; ++iU) - { - if (*(do_aik++)) - { - for (std::size_t i_cell = 0; i_cell < n_cells; ++i_cell) - U_vector[uik_nkj->first + i_cell] = A_vector[*aik + i_cell]; - ++aik; - } - for (std::size_t ikj = 0; ikj < uik_nkj->second; ++ikj) - { - for (std::size_t i_cell = 0; i_cell < n_cells; ++i_cell) - U_vector[uik_nkj->first + i_cell] -= L_vector[lij_ujk->first + i_cell] * U_vector[lij_ujk->second + i_cell]; - ++lij_ujk; - } - ++uik_nkj; - } - // Lower triangular matrix - for (std::size_t i_cell = 0; i_cell < n_cells; ++i_cell) - L_vector[lki_nkj->first + i_cell] = 1.0; - ++lki_nkj; - for (std::size_t iL = 0; iL < inLU.first; ++iL) - { - if (*(do_aki++)) - { - for (std::size_t i_cell = 0; i_cell < n_cells; ++i_cell) - L_vector[lki_nkj->first + i_cell] = A_vector[*aki + i_cell]; - ++aki; - } - for (std::size_t ikj = 0; ikj < lki_nkj->second; ++ikj) - { - for (std::size_t i_cell = 0; i_cell < n_cells; ++i_cell) - L_vector[lki_nkj->first + i_cell] -= L_vector[lkj_uji->first + i_cell] * U_vector[lkj_uji->second + i_cell]; - ++lkj_uji; - } - for (std::size_t i_cell = 0; i_cell < n_cells; ++i_cell) - L_vector[lki_nkj->first + i_cell] /= U_vector[*uii + i_cell]; - ++lki_nkj; - ++uii; - } - } - } - } } // namespace micm + +#include "lu_decomposition.inl" diff --git a/include/micm/solver/lu_decomposition.inl b/include/micm/solver/lu_decomposition.inl new file mode 100644 index 000000000..91aeebaf2 --- /dev/null +++ b/include/micm/solver/lu_decomposition.inl @@ -0,0 +1,266 @@ +namespace micm +{ + + inline LuDecomposition::LuDecomposition() + { + } + + template + inline LuDecomposition::LuDecomposition(const SparseMatrix& matrix) + { + std::size_t n = matrix[0].size(); + auto LU = GetLUMatrices(matrix, T{}); + const auto& L_row_start = LU.first.RowStartVector(); + const auto& L_row_ids = LU.first.RowIdsVector(); + const auto& U_row_start = LU.second.RowStartVector(); + const auto& U_row_ids = LU.second.RowIdsVector(); + for (std::size_t i = 0; i < matrix[0].size(); ++i) + { + std::pair iLU(0, 0); + // Upper triangular matrix + for (std::size_t k = i; k < n; ++k) + { + std::size_t nkj = 0; + for (std::size_t j_id = L_row_start[i]; j_id < L_row_start[i + 1]; ++j_id) + { + std::size_t j = L_row_ids[j_id]; + if (j >= i) + break; + if (LU.second.IsZero(j, k)) + continue; + ++nkj; + lij_ujk_.push_back(std::make_pair(LU.first.VectorIndex(0, i, j), LU.second.VectorIndex(0, j, k))); + } + if (matrix.IsZero(i, k)) + { + if (nkj == 0 && k != i) + continue; + do_aik_.push_back(false); + } + else + { + do_aik_.push_back(true); + aik_.push_back(matrix.VectorIndex(0, i, k)); + } + uik_nkj_.push_back(std::make_pair(LU.second.VectorIndex(0, i, k), nkj)); + ++(iLU.second); + } + // Lower triangular matrix + lki_nkj_.push_back(std::make_pair(LU.first.VectorIndex(0, i, i), 0)); + for (std::size_t k = i + 1; k < n; ++k) + { + std::size_t nkj = 0; + for (std::size_t j_id = L_row_start[k]; j_id < L_row_start[k + 1]; ++j_id) + { + std::size_t j = L_row_ids[j_id]; + if (j >= i) + break; + if (LU.second.IsZero(j, i)) + continue; + ++nkj; + lkj_uji_.push_back(std::make_pair(LU.first.VectorIndex(0, k, j), LU.second.VectorIndex(0, j, i))); + } + if (matrix.IsZero(k, i)) + { + if (nkj == 0) + continue; + do_aki_.push_back(false); + } + else + { + do_aki_.push_back(true); + aki_.push_back(matrix.VectorIndex(0, k, i)); + } + uii_.push_back(LU.second.VectorIndex(0, i, i)); + lki_nkj_.push_back(std::make_pair(LU.first.VectorIndex(0, k, i), nkj)); + ++(iLU.first); + } + niLU_.push_back(iLU); + } + } + + template + inline std::pair, SparseMatrix> LuDecomposition::GetLUMatrices( + const SparseMatrix& A, + T initial_value) + { + std::size_t n = A[0].size(); + std::set> L_ids, U_ids; + const auto& row_start = A.RowStartVector(); + const auto& row_ids = A.RowIdsVector(); + for (std::size_t i = 0; i < n; ++i) + { + // Upper triangular matrix + for (std::size_t k = i; k < n; ++k) + { + if (!A.IsZero(i, k) || k == i) + { + U_ids.insert(std::make_pair(i, k)); + continue; + } + for (std::size_t j = 0; j < i; ++j) + { + if (L_ids.find(std::make_pair(i, j)) != L_ids.end() && U_ids.find(std::make_pair(j, k)) != U_ids.end()) + { + U_ids.insert(std::make_pair(i, k)); + break; + } + } + } + // Lower triangular matrix + for (std::size_t k = i; k < n; ++k) + { + if (!A.IsZero(k, i) || k == i) + { + L_ids.insert(std::make_pair(k, i)); + continue; + } + for (std::size_t j = 0; j < i; ++j) + { + if (L_ids.find(std::make_pair(k, j)) != L_ids.end() && U_ids.find(std::make_pair(j, i)) != U_ids.end()) + { + L_ids.insert(std::make_pair(k, i)); + break; + } + } + } + } + auto L_builder = + micm::SparseMatrix::create(n).number_of_blocks(A.size()).initial_value(initial_value); + for (auto& pair : L_ids) + { + L_builder = L_builder.with_element(pair.first, pair.second); + } + auto U_builder = + micm::SparseMatrix::create(n).number_of_blocks(A.size()).initial_value(initial_value); + for (auto& pair : U_ids) + { + U_builder = U_builder.with_element(pair.first, pair.second); + } + std::pair, SparseMatrix> LU(L_builder, U_builder); + return LU; + } + + template class SparseMatrixPolicy> + requires(!VectorizableSparse>) inline void LuDecomposition::Decompose( + const SparseMatrixPolicy& A, + SparseMatrixPolicy& L, + SparseMatrixPolicy& U) const + { + // Loop over blocks + for (std::size_t i_block = 0; i_block < A.size(); ++i_block) + { + auto A_vector = std::next(A.AsVector().begin(), i_block * A.FlatBlockSize()); + auto L_vector = std::next(L.AsVector().begin(), i_block * L.FlatBlockSize()); + auto U_vector = std::next(U.AsVector().begin(), i_block * U.FlatBlockSize()); + auto do_aik = do_aik_.begin(); + auto aik = aik_.begin(); + auto uik_nkj = uik_nkj_.begin(); + auto lij_ujk = lij_ujk_.begin(); + auto do_aki = do_aki_.begin(); + auto aki = aki_.begin(); + auto lki_nkj = lki_nkj_.begin(); + auto lkj_uji = lkj_uji_.begin(); + auto uii = uii_.begin(); + for (auto& inLU : niLU_) + { + // Upper trianglur matrix + for (std::size_t iU = 0; iU < inLU.second; ++iU) + { + if (*(do_aik++)) + U_vector[uik_nkj->first] = A_vector[*(aik++)]; + for (std::size_t ikj = 0; ikj < uik_nkj->second; ++ikj) + { + U_vector[uik_nkj->first] -= L_vector[lij_ujk->first] * U_vector[lij_ujk->second]; + ++lij_ujk; + } + ++uik_nkj; + } + // Lower triangular matrix + L_vector[(lki_nkj++)->first] = 1.0; + for (std::size_t iL = 0; iL < inLU.first; ++iL) + { + if (*(do_aki++)) + L_vector[lki_nkj->first] = A_vector[*(aki++)]; + for (std::size_t ikj = 0; ikj < lki_nkj->second; ++ikj) + { + L_vector[lki_nkj->first] -= L_vector[lkj_uji->first] * U_vector[lkj_uji->second]; + ++lkj_uji; + } + L_vector[lki_nkj->first] /= U_vector[*uii]; + ++lki_nkj; + ++uii; + } + } + } + } + + template class SparseMatrixPolicy> + requires(VectorizableSparse>) inline void LuDecomposition::Decompose( + const SparseMatrixPolicy& A, + SparseMatrixPolicy& L, + SparseMatrixPolicy& U) const + { + const std::size_t n_cells = A.GroupVectorSize(); + // Loop over groups of blocks + for (std::size_t i_group = 0; i_group < A.NumberOfGroups(A.size()); ++i_group) + { + auto A_vector = std::next(A.AsVector().begin(), i_group * A.GroupSize(A.FlatBlockSize())); + auto L_vector = std::next(L.AsVector().begin(), i_group * L.GroupSize(L.FlatBlockSize())); + auto U_vector = std::next(U.AsVector().begin(), i_group * U.GroupSize(U.FlatBlockSize())); + auto do_aik = do_aik_.begin(); + auto aik = aik_.begin(); + auto uik_nkj = uik_nkj_.begin(); + auto lij_ujk = lij_ujk_.begin(); + auto do_aki = do_aki_.begin(); + auto aki = aki_.begin(); + auto lki_nkj = lki_nkj_.begin(); + auto lkj_uji = lkj_uji_.begin(); + auto uii = uii_.begin(); + for (auto& inLU : niLU_) + { + // Upper trianglur matrix + for (std::size_t iU = 0; iU < inLU.second; ++iU) + { + if (*(do_aik++)) + { + for (std::size_t i_cell = 0; i_cell < n_cells; ++i_cell) + U_vector[uik_nkj->first + i_cell] = A_vector[*aik + i_cell]; + ++aik; + } + for (std::size_t ikj = 0; ikj < uik_nkj->second; ++ikj) + { + for (std::size_t i_cell = 0; i_cell < n_cells; ++i_cell) + U_vector[uik_nkj->first + i_cell] -= L_vector[lij_ujk->first + i_cell] * U_vector[lij_ujk->second + i_cell]; + ++lij_ujk; + } + ++uik_nkj; + } + // Lower triangular matrix + for (std::size_t i_cell = 0; i_cell < n_cells; ++i_cell) + L_vector[lki_nkj->first + i_cell] = 1.0; + ++lki_nkj; + for (std::size_t iL = 0; iL < inLU.first; ++iL) + { + if (*(do_aki++)) + { + for (std::size_t i_cell = 0; i_cell < n_cells; ++i_cell) + L_vector[lki_nkj->first + i_cell] = A_vector[*aki + i_cell]; + ++aki; + } + for (std::size_t ikj = 0; ikj < lki_nkj->second; ++ikj) + { + for (std::size_t i_cell = 0; i_cell < n_cells; ++i_cell) + L_vector[lki_nkj->first + i_cell] -= L_vector[lkj_uji->first + i_cell] * U_vector[lkj_uji->second + i_cell]; + ++lkj_uji; + } + for (std::size_t i_cell = 0; i_cell < n_cells; ++i_cell) + L_vector[lki_nkj->first + i_cell] /= U_vector[*uii + i_cell]; + ++lki_nkj; + ++uii; + } + } + } + } + +} // namespace micm \ No newline at end of file diff --git a/include/micm/solver/rosenbrock.hpp b/include/micm/solver/rosenbrock.hpp index 2bd25ab89..c662fe456 100644 --- a/include/micm/solver/rosenbrock.hpp +++ b/include/micm/solver/rosenbrock.hpp @@ -82,52 +82,8 @@ namespace micm size_t number_of_grid_cells_{ 1 }; // Number of grid cells to solve simultaneously bool reorder_state_{ true }; // Reorder state during solver construction to minimize LU fill-in - void print() const - { - std::cout << "stages_: " << stages_ << std::endl; - std::cout << "upper_limit_tolerance_: " << upper_limit_tolerance_ << std::endl; - std::cout << "max_number_of_steps_: " << max_number_of_steps_ << std::endl; - std::cout << "round_off_: " << round_off_ << std::endl; - std::cout << "factor_min_: " << factor_min_ << std::endl; - std::cout << "factor_max_: " << factor_max_ << std::endl; - std::cout << "rejection_factor_decrease_: " << rejection_factor_decrease_ << std::endl; - std::cout << "safety_factor_: " << safety_factor_ << std::endl; - std::cout << "h_min_: " << h_min_ << std::endl; - std::cout << "h_max_: " << h_max_ << std::endl; - std::cout << "h_start_: " << h_start_ << std::endl; - std::cout << "new_function_evaluation_: "; - for (bool val : new_function_evaluation_) - std::cout << val << " "; - std::cout << std::endl; - std::cout << "estimator_of_local_order_: " << estimator_of_local_order_ << std::endl; - std::cout << "a_: "; - for (double val : a_) - std::cout << val << " "; - std::cout << std::endl; - std::cout << "c_: "; - for (double val : c_) - std::cout << val << " "; - std::cout << std::endl; - std::cout << "m_: "; - for (double val : m_) - std::cout << val << " "; - std::cout << std::endl; - std::cout << "e_: "; - for (double val : e_) - std::cout << val << " "; - std::cout << std::endl; - std::cout << "alpha_: "; - for (double val : alpha_) - std::cout << val << " "; - std::cout << std::endl; - std::cout << "gamma_: "; - for (double val : gamma_) - std::cout << val << " "; - std::cout << std::endl; - std::cout << "absolute_tolerance_: " << absolute_tolerance_ << std::endl; - std::cout << "relative_tolerance_: " << relative_tolerance_ << std::endl; - std::cout << "number_of_grid_cells_: " << number_of_grid_cells_ << std::endl; - } + // Print RosenbrockSolverParameters to console + void print() const; static RosenbrockSolverParameters two_stage_rosenbrock_parameters( size_t number_of_grid_cells = 1, @@ -150,316 +106,6 @@ namespace micm RosenbrockSolverParameters() = default; }; - RosenbrockSolverParameters RosenbrockSolverParameters::two_stage_rosenbrock_parameters( - size_t number_of_grid_cells, - bool reorder_state) - { - // an L-stable method, 2 stages, order 2 - - RosenbrockSolverParameters parameters; - double g = 1.0 + 1.0 / std::sqrt(2.0); - - parameters.stages_ = 2; - - parameters.a_.fill(0); - parameters.a_[0] = 1.0 / g; - - parameters.c_.fill(0); - parameters.c_[0] = -2.0 / g; - - // Both stages require a new function evaluation - parameters.new_function_evaluation_.fill(true); - - parameters.m_.fill(0); - parameters.m_[0] = (3.0) / (2.0 * g); - parameters.m_[1] = (1.0) / (2.0 * g); - - parameters.e_.fill(0); - parameters.e_[0] = 1.0 / (2.0 * g); - parameters.e_[1] = 1.0 / (2.0 * g); - - parameters.estimator_of_local_order_ = 2.0; - - parameters.alpha_.fill(0); - parameters.alpha_[0] = 0.0; - parameters.alpha_[1] = 1.0; - - parameters.gamma_.fill(0); - parameters.gamma_[0] = g; - parameters.gamma_[1] = -g; - - parameters.number_of_grid_cells_ = number_of_grid_cells; - parameters.reorder_state_ = reorder_state; - - return parameters; - } - - RosenbrockSolverParameters RosenbrockSolverParameters::three_stage_rosenbrock_parameters( - size_t number_of_grid_cells, - bool reorder_state) - { - // an L-stable method, 3 stages, order 3, 2 function evaluations - // - // original formaulation for three stages: - // Sandu, A., Verwer, J.G., Blom, J.G., Spee, E.J., Carmichael, G.R., Potra, F.A., 1997. - // Benchmarking stiff ode solvers for atmospheric chemistry problems II: Rosenbrock solvers. - // Atmospheric Environment 31, 3459–3472. https://doi.org/10.1016/S1352-2310(97)83212-8 - RosenbrockSolverParameters parameters; - - parameters.stages_ = 3; - - parameters.a_.fill(0); - parameters.a_[0] = 1; - parameters.a_[1] = 1; - parameters.a_[2] = 0; - - parameters.c_.fill(0); - parameters.c_[0] = -0.10156171083877702091975600115545e+01; - parameters.c_[1] = 0.40759956452537699824805835358067e+01; - parameters.c_[2] = 0.92076794298330791242156818474003e+01; - - parameters.new_function_evaluation_.fill(false); - parameters.new_function_evaluation_[0] = true; - parameters.new_function_evaluation_[1] = true; - parameters.new_function_evaluation_[2] = false; - parameters.m_.fill(0); - parameters.m_[0] = 0.1e+01; - parameters.m_[1] = 0.61697947043828245592553615689730e+01; - parameters.m_[2] = -0.42772256543218573326238373806514; - - parameters.e_.fill(0); - parameters.e_[0] = 0.5; - parameters.e_[1] = -0.29079558716805469821718236208017e+01; - parameters.e_[2] = 0.22354069897811569627360909276199; - - parameters.estimator_of_local_order_ = 3; - - parameters.alpha_.fill(0); - parameters.alpha_[0] = 0; - parameters.alpha_[1] = 0.43586652150845899941601945119356; - parameters.alpha_[2] = 0.43586652150845899941601945119356; - - parameters.gamma_.fill(0); - parameters.gamma_[0] = 0.43586652150845899941601945119356; - parameters.gamma_[1] = 0.24291996454816804366592249683314; - parameters.gamma_[2] = 0.21851380027664058511513169485832e+01; - - parameters.number_of_grid_cells_ = number_of_grid_cells; - parameters.reorder_state_ = reorder_state; - - return parameters; - } - - RosenbrockSolverParameters RosenbrockSolverParameters::four_stage_rosenbrock_parameters( - size_t number_of_grid_cells, - bool reorder_state) - { - // L-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 4 STAGES - // L-STABLE EMBEDDED ROSENBROCK METHOD OF ORDER 3 - // - // E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL - // EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS. - // SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS, - // SPRINGER-VERLAG (1990) - RosenbrockSolverParameters parameters; - - parameters.stages_ = 4; - - parameters.a_.fill(0); - parameters.a_[0] = 0.2000000000000000E+01; - parameters.a_[1] = 0.1867943637803922E+01; - parameters.a_[2] = 0.2344449711399156; - parameters.a_[3] = parameters.a_[1]; - parameters.a_[4] = parameters.a_[2]; - parameters.a_[5] = 0.0; - - parameters.c_.fill(0); - parameters.c_[0] = -0.7137615036412310E+01; - parameters.c_[1] = 0.2580708087951457E+01; - parameters.c_[2] = 0.6515950076447975; - parameters.c_[3] = -0.2137148994382534E+01; - parameters.c_[4] = -0.3214669691237626; - parameters.c_[5] = -0.6949742501781779; - - parameters.new_function_evaluation_.fill(false); - parameters.new_function_evaluation_[0] = true; - parameters.new_function_evaluation_[1] = true; - parameters.new_function_evaluation_[2] = true; - parameters.new_function_evaluation_[3] = false; - - parameters.m_.fill(0); - parameters.m_[0] = 0.2255570073418735E+01; - parameters.m_[1] = 0.2870493262186792; - parameters.m_[2] = 0.4353179431840180; - parameters.m_[3] = 0.1093502252409163E+01; - - parameters.e_.fill(0); - parameters.e_[0] = -0.2815431932141155; - parameters.e_[1] = -0.7276199124938920E-01; - parameters.e_[2] = -0.1082196201495311; - parameters.e_[3] = -0.1093502252409163E+01; - - parameters.estimator_of_local_order_ = 4.0; - - parameters.alpha_.fill(0); - parameters.alpha_[0] = 0.0; - parameters.alpha_[1] = 0.1145640000000000E+01; - parameters.alpha_[2] = 0.6552168638155900; - parameters.alpha_[3] = parameters.alpha_[2]; - - parameters.gamma_.fill(0); - parameters.gamma_[0] = 0.5728200000000000; - parameters.gamma_[1] = -0.1769193891319233E+01; - parameters.gamma_[2] = 0.7592633437920482; - parameters.gamma_[3] = -0.1049021087100450; - - parameters.number_of_grid_cells_ = number_of_grid_cells; - parameters.reorder_state_ = reorder_state; - - return parameters; - } - - RosenbrockSolverParameters RosenbrockSolverParameters::four_stage_differential_algebraic_rosenbrock_parameters( - size_t number_of_grid_cells, - bool reorder_state) - { - // A STIFFLY-STABLE METHOD, 4 stages, order 3 - RosenbrockSolverParameters parameters; - - // Set the number of stages - parameters.stages_ = 4; - - parameters.a_.fill(0.0); - parameters.a_[0] = 0.0; - parameters.a_[1] = 2.0; - parameters.a_[2] = 0.0; - parameters.a_[3] = 2.0; - parameters.a_[4] = 0.0; - parameters.a_[5] = 1.0; - - parameters.c_.fill(0.0); - parameters.c_[0] = 4.0; - parameters.c_[1] = 1.0; - parameters.c_[2] = -1.0; - parameters.c_[3] = 1.0; - parameters.c_[4] = -1.0; - parameters.c_[5] = -(8.0 / 3.0); - - parameters.new_function_evaluation_.fill(false); - parameters.new_function_evaluation_[0] = true; - parameters.new_function_evaluation_[2] = true; - parameters.new_function_evaluation_[3] = true; - - parameters.m_.fill(0.0); - parameters.m_[0] = 2.0; - parameters.m_[2] = 1.0; - parameters.m_[3] = 1.0; - - parameters.e_.fill(0.0); - parameters.e_[3] = 1.0; - - parameters.estimator_of_local_order_ = 3.0; - - parameters.alpha_.fill(0.0); - parameters.alpha_[2] = 1.0; - parameters.alpha_[3] = 1.0; - - parameters.gamma_.fill(0.0); - parameters.gamma_[0] = 0.5; - parameters.gamma_[1] = 1.5; - - parameters.number_of_grid_cells_ = number_of_grid_cells; - parameters.reorder_state_ = reorder_state; - - return parameters; - } - - RosenbrockSolverParameters RosenbrockSolverParameters::six_stage_differential_algebraic_rosenbrock_parameters( - size_t number_of_grid_cells, - bool reorder_state) - { - // STIFFLY-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 6 STAGES - // - // E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL - // EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS. - // SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS, - // SPRINGER-VERLAG (1996) - - RosenbrockSolverParameters parameters; - - parameters.stages_ = 6; - - parameters.alpha_.fill(0.0); - parameters.alpha_[0] = 0.000; - parameters.alpha_[1] = 0.386; - parameters.alpha_[2] = 0.210; - parameters.alpha_[3] = 0.630; - parameters.alpha_[4] = 1.000; - parameters.alpha_[5] = 1.000; - - parameters.gamma_.fill(0.0); - parameters.gamma_[0] = 0.2500000000000000; - parameters.gamma_[1] = -0.1043000000000000; - parameters.gamma_[2] = 0.1035000000000000; - parameters.gamma_[3] = -0.3620000000000023E-01; - parameters.gamma_[4] = 0.0; - parameters.gamma_[5] = 0.0; - - parameters.a_.fill(0.0); - parameters.a_[0] = 0.1544000000000000E+01; - parameters.a_[1] = 0.9466785280815826; - parameters.a_[2] = 0.2557011698983284; - parameters.a_[3] = 0.3314825187068521E+01; - parameters.a_[4] = 0.2896124015972201E+01; - parameters.a_[5] = 0.9986419139977817; - parameters.a_[6] = 0.1221224509226641E+01; - parameters.a_[7] = 0.6019134481288629E+01; - parameters.a_[8] = 0.1253708332932087E+02; - parameters.a_[9] = -0.6878860361058950; - parameters.a_[10] = parameters.a_[6]; - parameters.a_[11] = parameters.a_[7]; - parameters.a_[12] = parameters.a_[8]; - parameters.a_[13] = parameters.a_[9]; - parameters.a_[14] = 1.0; - - parameters.c_.fill(0.0); - parameters.c_[0] = -0.5668800000000000E+01; - parameters.c_[1] = -0.2430093356833875E+01; - parameters.c_[2] = -0.2063599157091915; - parameters.c_[3] = -0.1073529058151375; - parameters.c_[4] = -0.9594562251023355E+01; - parameters.c_[5] = -0.2047028614809616E+02; - parameters.c_[6] = 0.7496443313967647E+01; - parameters.c_[7] = -0.1024680431464352E+02; - parameters.c_[8] = -0.3399990352819905E+02; - parameters.c_[9] = 0.1170890893206160E+02; - parameters.c_[10] = 0.8083246795921522E+01; - parameters.c_[11] = -0.7981132988064893E+01; - parameters.c_[12] = -0.3152159432874371E+02; - parameters.c_[13] = 0.1631930543123136E+02; - parameters.c_[14] = -0.6058818238834054E+01; - - parameters.m_.fill(0.0); - parameters.m_[0] = parameters.a_[6]; - parameters.m_[1] = parameters.a_[7]; - parameters.m_[2] = parameters.a_[8]; - parameters.m_[3] = parameters.a_[9]; - parameters.m_[4] = 1.0; - parameters.m_[5] = 1.0; - - parameters.e_.fill(0.0); - parameters.e_[5] = 1.0; - - parameters.new_function_evaluation_.fill(true); - - parameters.estimator_of_local_order_ = 4.0; - - parameters.number_of_grid_cells_ = number_of_grid_cells; - parameters.reorder_state_ = reorder_state; - - return parameters; - } - /// @brief An implementation of the Chapman mechnanism solver /// /// The template parameter is the type of matrix to use @@ -553,10 +199,10 @@ namespace micm /// @brief compute [alpha * I - dforce_dy] /// @param jacobian Jacobian matrix (dforce_dy) /// @param alpha - inline void AlphaMinusJacobian(SparseMatrixPolicy& jacobian, const double& alpha) const - requires(!VectorizableSparse>); - inline void AlphaMinusJacobian(SparseMatrixPolicy& jacobian, const double& alpha) const - requires(VectorizableSparse>); + void AlphaMinusJacobian(SparseMatrixPolicy& jacobian, const double& alpha) const + requires(!VectorizableSparse>); + void AlphaMinusJacobian(SparseMatrixPolicy& jacobian, const double& alpha) const + requires(VectorizableSparse>); /// @brief Update the rate constants for the environment state /// @param state The current state of the chemical system @@ -566,7 +212,6 @@ namespace micm /// @param rate_constants List of rate constants for each needed species /// @param number_densities The number density of each species /// @param jacobian The matrix of partial derivatives - /// @return The jacobian void CalculateJacobian( const MatrixPolicy& rate_constants, const MatrixPolicy& number_densities, @@ -594,403 +239,7 @@ namespace micm double NormalizedError(const MatrixPolicy& y, const MatrixPolicy& y_new, const MatrixPolicy& errors); }; + +} // namespace micm - template class MatrixPolicy, template class SparseMatrixPolicy> - void RosenbrockSolver::SolverStats::Reset() - { - function_calls = 0; - jacobian_updates = 0; - number_of_steps = 0; - accepted = 0; - rejected = 0; - decompositions = 0; - solves = 0; - singular = 0; - total_steps = 0; - } - - template class MatrixPolicy, template class SparseMatrixPolicy> - std::string RosenbrockSolver::SolverStats::State( - const RosenbrockSolver::SolverState& state) const - { - switch (state) - { - case RosenbrockSolver::SolverState::NotYetCalled: return "Not Yet Called"; - case RosenbrockSolver::SolverState::Running: return "Running"; - case RosenbrockSolver::SolverState::Converged: return "Converged"; - case RosenbrockSolver::SolverState::ConvergenceExceededMaxSteps: - return "Convergence Exceeded Max Steps"; - case RosenbrockSolver::SolverState::StepSizeTooSmall: return "Step Size Too Small"; - case RosenbrockSolver::SolverState::RepeatedlySingularMatrix: - return "Repeatedly Singular Matrix"; - case RosenbrockSolver::SolverState::NaNDetected: return "NaNDetected"; - default: return "Unknown"; - } - return ""; - } - - template class MatrixPolicy, template class SparseMatrixPolicy> - inline RosenbrockSolver::RosenbrockSolver() - : system_(), - processes_(), - parameters_(RosenbrockSolverParameters::three_stage_rosenbrock_parameters()), - process_set_(), - stats_(), - jacobian_(), - linear_solver_(), - jacobian_diagonal_elements_(), - N_(system_.StateSize() * parameters_.number_of_grid_cells_) - { - } - - template class MatrixPolicy, template class SparseMatrixPolicy> - inline RosenbrockSolver::RosenbrockSolver( - const System& system, - const std::vector& processes, - const RosenbrockSolverParameters& parameters) - : system_(system), - processes_(processes), - parameters_(parameters), - state_reordering_(), - process_set_(), - stats_(), - jacobian_(), - linear_solver_(), - jacobian_diagonal_elements_(), - N_(system_.StateSize() * parameters_.number_of_grid_cells_) - { - // generate a state-vector reordering function to reduce fill-in in linear solver - if (parameters_.reorder_state_) - { - // get unsorted Jacobian non-zero elements - auto unsorted_process_set = ProcessSet(processes, GetState()); - auto unsorted_jac_elements = unsorted_process_set.NonZeroJacobianElements(); - MatrixPolicy unsorted_jac_non_zeros(system_.StateSize(), system_.StateSize(), 0); - for (auto& elem : unsorted_jac_elements) - unsorted_jac_non_zeros[elem.first][elem.second] = 1; - auto reorder_map = DiagonalMarkowitzReorder(unsorted_jac_non_zeros); - state_reordering_ = [=](const std::vector& variables, const std::size_t i) - { return variables[reorder_map[i]]; }; - } - process_set_ = ProcessSet(processes, GetState()); - auto builder = - SparseMatrixPolicy::create(system_.StateSize()).number_of_blocks(parameters_.number_of_grid_cells_); - auto jac_elements = process_set_.NonZeroJacobianElements(); - for (auto& elem : jac_elements) - builder = builder.with_element(elem.first, elem.second); - // Always include diagonal elements - for (std::size_t i = 0; i < system_.StateSize(); ++i) - builder = builder.with_element(i, i); - - jacobian_ = builder; - linear_solver_ = LinearSolver(jacobian_, 1.0e-30); - process_set_.SetJacobianFlatIds(jacobian_); - for (std::size_t i = 0; i < jacobian_[0].size(); ++i) - jacobian_diagonal_elements_.push_back(jacobian_.VectorIndex(0, i, i)); - } - - template class MatrixPolicy, template class SparseMatrixPolicy> - inline RosenbrockSolver::~RosenbrockSolver() - { - } - - template class MatrixPolicy, template class SparseMatrixPolicy> - inline State RosenbrockSolver::GetState() const - { - std::vector param_labels{}; - for (const auto& process : processes_) - if (process.rate_constant_) - for (auto& label : process.rate_constant_->CustomParameters()) - param_labels.push_back(label); - return State{ micm::StateParameters{ .state_variable_names_ = system_.UniqueNames(state_reordering_), - .custom_rate_parameter_labels_ = param_labels, - .number_of_grid_cells_ = parameters_.number_of_grid_cells_, - .number_of_rate_constants_ = processes_.size() } }; - } - - template class MatrixPolicy, template class SparseMatrixPolicy> - inline typename RosenbrockSolver::SolverResult - RosenbrockSolver::Solve(double time_step, State& state) noexcept - { - typename RosenbrockSolver::SolverResult result{}; - result.state_ = SolverState::Running; - MatrixPolicy Y(state.variables_); - MatrixPolicy Ynew(Y.size(), Y[0].size(), 0.0); - MatrixPolicy initial_forcing(Y.size(), Y[0].size(), 0.0); - MatrixPolicy forcing(Y.size(), Y[0].size(), 0.0); - MatrixPolicy temp(Y.size(), Y[0].size(), 0.0); - std::vector> K{}; - - parameters_.h_max_ = time_step; - parameters_.h_start_ = std::max(parameters_.h_min_, delta_min_); - - stats_.Reset(); - UpdateState(state); - - for (std::size_t i = 0; i < parameters_.stages_; ++i) - K.push_back(MatrixPolicy(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_)); - - if (std::abs(H) <= 10 * parameters_.round_off_) - H = delta_min_; - - bool reject_last_h = false; - bool reject_more_h = false; - - while ((present_time - time_step + parameters_.round_off_) <= 0 && (result.state_ == SolverState::Running)) - { - if (stats_.number_of_steps > parameters_.max_number_of_steps_) - { - result.state_ = SolverState::ConvergenceExceededMaxSteps; - break; - } - - if (((present_time + 0.1 * H) == present_time) || (H <= parameters_.round_off_)) - { - result.state_ = SolverState::StepSizeTooSmall; - break; - } - - // Limit H if necessary to avoid going beyond the specified chemistry time step - H = std::min(H, std::abs(time_step - present_time)); - - // compute the concentrations at the current time - CalculateForcing(state.rate_constants_, Y, initial_forcing); - - // compute the jacobian at the current time - CalculateJacobian(state.rate_constants_, Y, jacobian_); - - bool accepted = false; - // Repeat step calculation until current step accepted - while (!accepted) - { - bool is_singular{ false }; - // Form and factor the rosenbrock ode jacobian - LinearFactor(H, parameters_.gamma_[0], is_singular, Y, state.rate_constants_); - if (is_singular) - { - result.state_ = SolverState::RepeatedlySingularMatrix; - break; - } - - // Compute the stages - for (uint64_t stage = 0; stage < parameters_.stages_; ++stage) - { - double stage_combinations = ((stage + 1) - 1) * ((stage + 1) - 2) / 2; - if (stage == 0) - { - forcing = initial_forcing; - } - else - { - if (parameters_.new_function_evaluation_[stage]) - { - Ynew.AsVector().assign(Y.AsVector().begin(), Y.AsVector().end()); - for (uint64_t j = 0; j < stage; ++j) - { - auto a = parameters_.a_[stage_combinations + j]; - Ynew.ForEach([&](double& iYnew, double& iKj) { iYnew += a * iKj; }, K[j]); - } - CalculateForcing(state.rate_constants_, Ynew, forcing); - } - } - K[stage].AsVector().assign(forcing.AsVector().begin(), forcing.AsVector().end()); - for (uint64_t j = 0; j < stage; ++j) - { - auto HC = parameters_.c_[stage_combinations + j] / H; - K[stage].ForEach([&](double& iKstage, double& iKj) { iKstage += HC * iKj; }, K[j]); - } - temp.AsVector().assign(K[stage].AsVector().begin(), K[stage].AsVector().end()); - linear_solver_.template Solve(temp, K[stage]); - stats_.solves += 1; - } - - // Compute the new solution - Ynew.AsVector().assign(Y.AsVector().begin(), Y.AsVector().end()); - for (uint64_t stage = 0; stage < parameters_.stages_; ++stage) - Ynew.ForEach([&](double& iYnew, double& iKstage) { iYnew += parameters_.m_[stage] * iKstage; }, K[stage]); - - // Compute the error estimation - MatrixPolicy Yerror(Y.size(), Y[0].size(), 0); - for (uint64_t stage = 0; stage < parameters_.stages_; ++stage) - Yerror.ForEach([&](double& iYerror, double& iKstage) { iYerror += parameters_.e_[stage] * iKstage; }, K[stage]); - - auto error = NormalizedError(Y, Ynew, Yerror); - - // New step size is bounded by FacMin <= Hnew/H <= FacMax - double fac = std::min( - parameters_.factor_max_, - std::max( - parameters_.factor_min_, - parameters_.safety_factor_ / std::pow(error, 1 / parameters_.estimator_of_local_order_))); - double Hnew = H * fac; - - // Check the error magnitude and adjust step size - stats_.number_of_steps += 1; - stats_.total_steps += 1; - - if (std::isnan(error)) - { - Y.AsVector().assign(Ynew.AsVector().begin(), Ynew.AsVector().end()); - result.state_ = SolverState::NaNDetected; - break; - } - else if ((error < 1) || (H < parameters_.h_min_)) - { - 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_)); - if (reject_last_h) - { - // No step size increase after a rejected step - Hnew = std::min(Hnew, H); - } - reject_last_h = false; - reject_more_h = false; - H = Hnew; - accepted = true; - } - else - { - // Reject step - if (reject_more_h) - { - Hnew = H * parameters_.rejection_factor_decrease_; - } - reject_more_h = reject_last_h; - reject_last_h = true; - H = Hnew; - if (stats_.accepted >= 1) - { - stats_.rejected += 1; - } - } - } - } - - result.state_ = SolverState::Converged; - result.final_time_ = present_time; - result.stats_ = stats_; - result.result_ = std::move(Y); - return result; - } - - template class MatrixPolicy, template class SparseMatrixPolicy> - inline void RosenbrockSolver::CalculateForcing( - const MatrixPolicy& rate_constants, - const MatrixPolicy& number_densities, - MatrixPolicy& forcing) - { - std::fill(forcing.AsVector().begin(), forcing.AsVector().end(), 0.0); - process_set_.AddForcingTerms(rate_constants, number_densities, forcing); - stats_.function_calls += 1; - } - - template class MatrixPolicy, template class SparseMatrixPolicy> - inline void RosenbrockSolver::AlphaMinusJacobian( - SparseMatrixPolicy& jacobian, - const double& alpha) const requires(!VectorizableSparse>) - { - for (auto& elem : jacobian.AsVector()) - elem = -elem; - for (std::size_t i_block = 0; i_block < jacobian.size(); ++i_block) - { - auto jacobian_vector = std::next(jacobian.AsVector().begin(), i_block * jacobian.FlatBlockSize()); - for (const auto& i_elem : jacobian_diagonal_elements_) - jacobian_vector[i_elem] += alpha; - } - } - - template class MatrixPolicy, template class SparseMatrixPolicy> - inline void RosenbrockSolver::AlphaMinusJacobian( - SparseMatrixPolicy& jacobian, - const double& alpha) const requires(VectorizableSparse>) - { - const std::size_t n_cells = jacobian.GroupVectorSize(); - for (auto& elem : jacobian.AsVector()) - elem = -elem; - for (std::size_t i_group = 0; i_group < jacobian.NumberOfGroups(jacobian.size()); ++i_group) - { - auto jacobian_vector = std::next(jacobian.AsVector().begin(), i_group * jacobian.GroupSize(jacobian.FlatBlockSize())); - for (const auto& i_elem : jacobian_diagonal_elements_) - for (std::size_t i_cell = 0; i_cell < n_cells; ++i_cell) - jacobian_vector[i_elem + i_cell] += alpha; - } - } - - template class MatrixPolicy, template class SparseMatrixPolicy> - inline void RosenbrockSolver::CalculateJacobian( - const MatrixPolicy& rate_constants, - const MatrixPolicy& number_densities, - SparseMatrixPolicy& jacobian) - { - std::fill(jacobian.AsVector().begin(), jacobian.AsVector().end(), 0.0); - process_set_.AddJacobianTerms(rate_constants, number_densities, jacobian); - stats_.jacobian_updates += 1; - } - - template class MatrixPolicy, template class SparseMatrixPolicy> - inline void RosenbrockSolver::UpdateState(State& state) - { - Process::UpdateState(processes_, state); - } - - template class MatrixPolicy, template class SparseMatrixPolicy> - inline void RosenbrockSolver::LinearFactor( - double& H, - const double gamma, - bool& singular, - const MatrixPolicy& number_densities, - const MatrixPolicy& rate_constants) - { - // TODO: invesitage this function. The fortran equivalent appears to have a bug. - // From my understanding the fortran do loop would only ever do one iteration and is equivalent to what's below - SparseMatrixPolicy jacobian = jacobian_; - uint64_t n_consecutive = 0; - singular = true; - while (true) - { - double alpha = 1 / (H * gamma); - AlphaMinusJacobian(jacobian, alpha); - linear_solver_.Factor(jacobian); - singular = false; // TODO This should be evaluated in some way - stats_.decompositions += 1; - if (!singular) - break; - stats_.singular += 1; - if (++n_consecutive > 5) - break; - H /= 2; - } - } - - template class MatrixPolicy, template class SparseMatrixPolicy> - inline double RosenbrockSolver::NormalizedError( - const MatrixPolicy& Y, - const MatrixPolicy& Ynew, - const MatrixPolicy& errors) - { - // Solving Ordinary Differential Equations II, page 123 - // https://link-springer-com.cuucar.idm.oclc.org/book/10.1007/978-3-642-05221-7 - - auto _y = Y.AsVector(); - auto _ynew = Ynew.AsVector(); - auto _errors = errors.AsVector(); - - double error = 0; - - for (size_t i = 0; i < N_; ++i) - { - double ymax = std::max(std::abs(_y[i]), std::abs(_ynew[i])); - double scale = parameters_.absolute_tolerance_ + parameters_.relative_tolerance_ * ymax; - error += std::pow(_errors[i] / scale, 2); - } - - double error_min_ = 1.0e-10; - return std::max(std::sqrt(error / N_), error_min_); - } -} // namespace micm \ No newline at end of file +#include "rosenbrock.inl" \ No newline at end of file diff --git a/include/micm/solver/rosenbrock.inl b/include/micm/solver/rosenbrock.inl new file mode 100644 index 000000000..64aa48d30 --- /dev/null +++ b/include/micm/solver/rosenbrock.inl @@ -0,0 +1,768 @@ +namespace micm +{ + // + // RosenbrockSolverParameters + // + void RosenbrockSolverParameters::print() const + { + std::cout << "stages_: " << stages_ << std::endl; + std::cout << "upper_limit_tolerance_: " << upper_limit_tolerance_ << std::endl; + std::cout << "max_number_of_steps_: " << max_number_of_steps_ << std::endl; + std::cout << "round_off_: " << round_off_ << std::endl; + std::cout << "factor_min_: " << factor_min_ << std::endl; + std::cout << "factor_max_: " << factor_max_ << std::endl; + std::cout << "rejection_factor_decrease_: " << rejection_factor_decrease_ << std::endl; + std::cout << "safety_factor_: " << safety_factor_ << std::endl; + std::cout << "h_min_: " << h_min_ << std::endl; + std::cout << "h_max_: " << h_max_ << std::endl; + std::cout << "h_start_: " << h_start_ << std::endl; + std::cout << "new_function_evaluation_: "; + for (bool val : new_function_evaluation_) + std::cout << val << " "; + std::cout << std::endl; + std::cout << "estimator_of_local_order_: " << estimator_of_local_order_ << std::endl; + std::cout << "a_: "; + for (double val : a_) + std::cout << val << " "; + std::cout << std::endl; + std::cout << "c_: "; + for (double val : c_) + std::cout << val << " "; + std::cout << std::endl; + std::cout << "m_: "; + for (double val : m_) + std::cout << val << " "; + std::cout << std::endl; + std::cout << "e_: "; + for (double val : e_) + std::cout << val << " "; + std::cout << std::endl; + std::cout << "alpha_: "; + for (double val : alpha_) + std::cout << val << " "; + std::cout << std::endl; + std::cout << "gamma_: "; + for (double val : gamma_) + std::cout << val << " "; + std::cout << std::endl; + std::cout << "absolute_tolerance_: " << absolute_tolerance_ << std::endl; + std::cout << "relative_tolerance_: " << relative_tolerance_ << std::endl; + std::cout << "number_of_grid_cells_: " << number_of_grid_cells_ << std::endl; + } + + RosenbrockSolverParameters RosenbrockSolverParameters::two_stage_rosenbrock_parameters( + size_t number_of_grid_cells, + bool reorder_state) + { + // an L-stable method, 2 stages, order 2 + + RosenbrockSolverParameters parameters; + double g = 1.0 + 1.0 / std::sqrt(2.0); + + parameters.stages_ = 2; + + parameters.a_.fill(0); + parameters.a_[0] = 1.0 / g; + + parameters.c_.fill(0); + parameters.c_[0] = -2.0 / g; + + // Both stages require a new function evaluation + parameters.new_function_evaluation_.fill(true); + + parameters.m_.fill(0); + parameters.m_[0] = (3.0) / (2.0 * g); + parameters.m_[1] = (1.0) / (2.0 * g); + + parameters.e_.fill(0); + parameters.e_[0] = 1.0 / (2.0 * g); + parameters.e_[1] = 1.0 / (2.0 * g); + + parameters.estimator_of_local_order_ = 2.0; + + parameters.alpha_.fill(0); + parameters.alpha_[0] = 0.0; + parameters.alpha_[1] = 1.0; + + parameters.gamma_.fill(0); + parameters.gamma_[0] = g; + parameters.gamma_[1] = -g; + + parameters.number_of_grid_cells_ = number_of_grid_cells; + parameters.reorder_state_ = reorder_state; + + return parameters; + } + + RosenbrockSolverParameters RosenbrockSolverParameters::three_stage_rosenbrock_parameters( + size_t number_of_grid_cells, + bool reorder_state) + { + // an L-stable method, 3 stages, order 3, 2 function evaluations + // + // original formaulation for three stages: + // Sandu, A., Verwer, J.G., Blom, J.G., Spee, E.J., Carmichael, G.R., Potra, F.A., 1997. + // Benchmarking stiff ode solvers for atmospheric chemistry problems II: Rosenbrock solvers. + // Atmospheric Environment 31, 3459–3472. https://doi.org/10.1016/S1352-2310(97)83212-8 + RosenbrockSolverParameters parameters; + + parameters.stages_ = 3; + + parameters.a_.fill(0); + parameters.a_[0] = 1; + parameters.a_[1] = 1; + parameters.a_[2] = 0; + + parameters.c_.fill(0); + parameters.c_[0] = -0.10156171083877702091975600115545e+01; + parameters.c_[1] = 0.40759956452537699824805835358067e+01; + parameters.c_[2] = 0.92076794298330791242156818474003e+01; + + parameters.new_function_evaluation_.fill(false); + parameters.new_function_evaluation_[0] = true; + parameters.new_function_evaluation_[1] = true; + parameters.new_function_evaluation_[2] = false; + parameters.m_.fill(0); + parameters.m_[0] = 0.1e+01; + parameters.m_[1] = 0.61697947043828245592553615689730e+01; + parameters.m_[2] = -0.42772256543218573326238373806514; + + parameters.e_.fill(0); + parameters.e_[0] = 0.5; + parameters.e_[1] = -0.29079558716805469821718236208017e+01; + parameters.e_[2] = 0.22354069897811569627360909276199; + + parameters.estimator_of_local_order_ = 3; + + parameters.alpha_.fill(0); + parameters.alpha_[0] = 0; + parameters.alpha_[1] = 0.43586652150845899941601945119356; + parameters.alpha_[2] = 0.43586652150845899941601945119356; + + parameters.gamma_.fill(0); + parameters.gamma_[0] = 0.43586652150845899941601945119356; + parameters.gamma_[1] = 0.24291996454816804366592249683314; + parameters.gamma_[2] = 0.21851380027664058511513169485832e+01; + + parameters.number_of_grid_cells_ = number_of_grid_cells; + parameters.reorder_state_ = reorder_state; + + return parameters; + } + + RosenbrockSolverParameters RosenbrockSolverParameters::four_stage_rosenbrock_parameters( + size_t number_of_grid_cells, + bool reorder_state) + { + // L-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 4 STAGES + // L-STABLE EMBEDDED ROSENBROCK METHOD OF ORDER 3 + // + // E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL + // EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS. + // SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS, + // SPRINGER-VERLAG (1990) + RosenbrockSolverParameters parameters; + + parameters.stages_ = 4; + + parameters.a_.fill(0); + parameters.a_[0] = 0.2000000000000000E+01; + parameters.a_[1] = 0.1867943637803922E+01; + parameters.a_[2] = 0.2344449711399156; + parameters.a_[3] = parameters.a_[1]; + parameters.a_[4] = parameters.a_[2]; + parameters.a_[5] = 0.0; + + parameters.c_.fill(0); + parameters.c_[0] = -0.7137615036412310E+01; + parameters.c_[1] = 0.2580708087951457E+01; + parameters.c_[2] = 0.6515950076447975; + parameters.c_[3] = -0.2137148994382534E+01; + parameters.c_[4] = -0.3214669691237626; + parameters.c_[5] = -0.6949742501781779; + + parameters.new_function_evaluation_.fill(false); + parameters.new_function_evaluation_[0] = true; + parameters.new_function_evaluation_[1] = true; + parameters.new_function_evaluation_[2] = true; + parameters.new_function_evaluation_[3] = false; + + parameters.m_.fill(0); + parameters.m_[0] = 0.2255570073418735E+01; + parameters.m_[1] = 0.2870493262186792; + parameters.m_[2] = 0.4353179431840180; + parameters.m_[3] = 0.1093502252409163E+01; + + parameters.e_.fill(0); + parameters.e_[0] = -0.2815431932141155; + parameters.e_[1] = -0.7276199124938920E-01; + parameters.e_[2] = -0.1082196201495311; + parameters.e_[3] = -0.1093502252409163E+01; + + parameters.estimator_of_local_order_ = 4.0; + + parameters.alpha_.fill(0); + parameters.alpha_[0] = 0.0; + parameters.alpha_[1] = 0.1145640000000000E+01; + parameters.alpha_[2] = 0.6552168638155900; + parameters.alpha_[3] = parameters.alpha_[2]; + + parameters.gamma_.fill(0); + parameters.gamma_[0] = 0.5728200000000000; + parameters.gamma_[1] = -0.1769193891319233E+01; + parameters.gamma_[2] = 0.7592633437920482; + parameters.gamma_[3] = -0.1049021087100450; + + parameters.number_of_grid_cells_ = number_of_grid_cells; + parameters.reorder_state_ = reorder_state; + + return parameters; + } + + RosenbrockSolverParameters RosenbrockSolverParameters::four_stage_differential_algebraic_rosenbrock_parameters( + size_t number_of_grid_cells, + bool reorder_state) + { + // A STIFFLY-STABLE METHOD, 4 stages, order 3 + RosenbrockSolverParameters parameters; + + // Set the number of stages + parameters.stages_ = 4; + + parameters.a_.fill(0.0); + parameters.a_[0] = 0.0; + parameters.a_[1] = 2.0; + parameters.a_[2] = 0.0; + parameters.a_[3] = 2.0; + parameters.a_[4] = 0.0; + parameters.a_[5] = 1.0; + + parameters.c_.fill(0.0); + parameters.c_[0] = 4.0; + parameters.c_[1] = 1.0; + parameters.c_[2] = -1.0; + parameters.c_[3] = 1.0; + parameters.c_[4] = -1.0; + parameters.c_[5] = -(8.0 / 3.0); + + parameters.new_function_evaluation_.fill(false); + parameters.new_function_evaluation_[0] = true; + parameters.new_function_evaluation_[2] = true; + parameters.new_function_evaluation_[3] = true; + + parameters.m_.fill(0.0); + parameters.m_[0] = 2.0; + parameters.m_[2] = 1.0; + parameters.m_[3] = 1.0; + + parameters.e_.fill(0.0); + parameters.e_[3] = 1.0; + + parameters.estimator_of_local_order_ = 3.0; + + parameters.alpha_.fill(0.0); + parameters.alpha_[2] = 1.0; + parameters.alpha_[3] = 1.0; + + parameters.gamma_.fill(0.0); + parameters.gamma_[0] = 0.5; + parameters.gamma_[1] = 1.5; + + parameters.number_of_grid_cells_ = number_of_grid_cells; + parameters.reorder_state_ = reorder_state; + + return parameters; + } + + RosenbrockSolverParameters RosenbrockSolverParameters::six_stage_differential_algebraic_rosenbrock_parameters( + size_t number_of_grid_cells, + bool reorder_state) + { + // STIFFLY-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 6 STAGES + // + // E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL + // EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS. + // SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS, + // SPRINGER-VERLAG (1996) + + RosenbrockSolverParameters parameters; + + parameters.stages_ = 6; + + parameters.alpha_.fill(0.0); + parameters.alpha_[0] = 0.000; + parameters.alpha_[1] = 0.386; + parameters.alpha_[2] = 0.210; + parameters.alpha_[3] = 0.630; + parameters.alpha_[4] = 1.000; + parameters.alpha_[5] = 1.000; + + parameters.gamma_.fill(0.0); + parameters.gamma_[0] = 0.2500000000000000; + parameters.gamma_[1] = -0.1043000000000000; + parameters.gamma_[2] = 0.1035000000000000; + parameters.gamma_[3] = -0.3620000000000023E-01; + parameters.gamma_[4] = 0.0; + parameters.gamma_[5] = 0.0; + + parameters.a_.fill(0.0); + parameters.a_[0] = 0.1544000000000000E+01; + parameters.a_[1] = 0.9466785280815826; + parameters.a_[2] = 0.2557011698983284; + parameters.a_[3] = 0.3314825187068521E+01; + parameters.a_[4] = 0.2896124015972201E+01; + parameters.a_[5] = 0.9986419139977817; + parameters.a_[6] = 0.1221224509226641E+01; + parameters.a_[7] = 0.6019134481288629E+01; + parameters.a_[8] = 0.1253708332932087E+02; + parameters.a_[9] = -0.6878860361058950; + parameters.a_[10] = parameters.a_[6]; + parameters.a_[11] = parameters.a_[7]; + parameters.a_[12] = parameters.a_[8]; + parameters.a_[13] = parameters.a_[9]; + parameters.a_[14] = 1.0; + + parameters.c_.fill(0.0); + parameters.c_[0] = -0.5668800000000000E+01; + parameters.c_[1] = -0.2430093356833875E+01; + parameters.c_[2] = -0.2063599157091915; + parameters.c_[3] = -0.1073529058151375; + parameters.c_[4] = -0.9594562251023355E+01; + parameters.c_[5] = -0.2047028614809616E+02; + parameters.c_[6] = 0.7496443313967647E+01; + parameters.c_[7] = -0.1024680431464352E+02; + parameters.c_[8] = -0.3399990352819905E+02; + parameters.c_[9] = 0.1170890893206160E+02; + parameters.c_[10] = 0.8083246795921522E+01; + parameters.c_[11] = -0.7981132988064893E+01; + parameters.c_[12] = -0.3152159432874371E+02; + parameters.c_[13] = 0.1631930543123136E+02; + parameters.c_[14] = -0.6058818238834054E+01; + + parameters.m_.fill(0.0); + parameters.m_[0] = parameters.a_[6]; + parameters.m_[1] = parameters.a_[7]; + parameters.m_[2] = parameters.a_[8]; + parameters.m_[3] = parameters.a_[9]; + parameters.m_[4] = 1.0; + parameters.m_[5] = 1.0; + + parameters.e_.fill(0.0); + parameters.e_[5] = 1.0; + + parameters.new_function_evaluation_.fill(true); + + parameters.estimator_of_local_order_ = 4.0; + + parameters.number_of_grid_cells_ = number_of_grid_cells; + parameters.reorder_state_ = reorder_state; + + return parameters; + } + + // + // RosenbrockSolver + // + template class MatrixPolicy, template class SparseMatrixPolicy> + inline void RosenbrockSolver::SolverStats::Reset() + { + function_calls = 0; + jacobian_updates = 0; + number_of_steps = 0; + accepted = 0; + rejected = 0; + decompositions = 0; + solves = 0; + singular = 0; + total_steps = 0; + } + + template class MatrixPolicy, template class SparseMatrixPolicy> + inline std::string RosenbrockSolver::SolverStats::State( + const RosenbrockSolver::SolverState& state) const + { + switch (state) + { + case RosenbrockSolver::SolverState::NotYetCalled: return "Not Yet Called"; + case RosenbrockSolver::SolverState::Running: return "Running"; + case RosenbrockSolver::SolverState::Converged: return "Converged"; + case RosenbrockSolver::SolverState::ConvergenceExceededMaxSteps: + return "Convergence Exceeded Max Steps"; + case RosenbrockSolver::SolverState::StepSizeTooSmall: return "Step Size Too Small"; + case RosenbrockSolver::SolverState::RepeatedlySingularMatrix: + return "Repeatedly Singular Matrix"; + case RosenbrockSolver::SolverState::NaNDetected: return "NaNDetected"; + default: return "Unknown"; + } + return ""; + } + + template class MatrixPolicy, template class SparseMatrixPolicy> + inline RosenbrockSolver::RosenbrockSolver() + : system_(), + processes_(), + parameters_(RosenbrockSolverParameters::three_stage_rosenbrock_parameters()), + state_reordering_(), + process_set_(), + stats_(), + jacobian_(), + linear_solver_(), + jacobian_diagonal_elements_(), + N_(system_.StateSize() * parameters_.number_of_grid_cells_) + { + } + + template class MatrixPolicy, template class SparseMatrixPolicy> + inline RosenbrockSolver::RosenbrockSolver( + const System& system, + const std::vector& processes, + const RosenbrockSolverParameters& parameters) + : system_(system), + processes_(processes), + parameters_(parameters), + state_reordering_(), + process_set_(), + stats_(), + jacobian_(), + linear_solver_(), + jacobian_diagonal_elements_(), + N_(system_.StateSize() * parameters_.number_of_grid_cells_) + { + // generate a state-vector reordering function to reduce fill-in in linear solver + if (parameters_.reorder_state_) + { + // get unsorted Jacobian non-zero elements + auto unsorted_process_set = ProcessSet(processes, GetState()); + auto unsorted_jac_elements = unsorted_process_set.NonZeroJacobianElements(); + MatrixPolicy unsorted_jac_non_zeros(system_.StateSize(), system_.StateSize(), 0); + for (auto& elem : unsorted_jac_elements) + unsorted_jac_non_zeros[elem.first][elem.second] = 1; + auto reorder_map = DiagonalMarkowitzReorder(unsorted_jac_non_zeros); + state_reordering_ = [=](const std::vector& variables, const std::size_t i) + { return variables[reorder_map[i]]; }; + } + process_set_ = ProcessSet(processes, GetState()); + auto builder = + SparseMatrixPolicy::create(system_.StateSize()).number_of_blocks(parameters_.number_of_grid_cells_); + auto jac_elements = process_set_.NonZeroJacobianElements(); + for (auto& elem : jac_elements) + builder = builder.with_element(elem.first, elem.second); + // Always include diagonal elements + for (std::size_t i = 0; i < system_.StateSize(); ++i) + builder = builder.with_element(i, i); + + jacobian_ = builder; + linear_solver_ = LinearSolver(jacobian_, 1.0e-30); + process_set_.SetJacobianFlatIds(jacobian_); + for (std::size_t i = 0; i < jacobian_[0].size(); ++i) + jacobian_diagonal_elements_.push_back(jacobian_.VectorIndex(0, i, i)); + } + + template class MatrixPolicy, template class SparseMatrixPolicy> + inline RosenbrockSolver::~RosenbrockSolver() + { + } + + template class MatrixPolicy, template class SparseMatrixPolicy> + inline State RosenbrockSolver::GetState() const + { + std::vector param_labels{}; + for (const auto& process : processes_) + if (process.rate_constant_) + for (auto& label : process.rate_constant_->CustomParameters()) + param_labels.push_back(label); + return State{ micm::StateParameters{ .state_variable_names_ = system_.UniqueNames(state_reordering_), + .custom_rate_parameter_labels_ = param_labels, + .number_of_grid_cells_ = parameters_.number_of_grid_cells_, + .number_of_rate_constants_ = processes_.size() } }; + } + + template class MatrixPolicy, template class SparseMatrixPolicy> + inline typename RosenbrockSolver::SolverResult + RosenbrockSolver::Solve(double time_step, State& state) noexcept + { + typename RosenbrockSolver::SolverResult result{}; + result.state_ = SolverState::Running; + MatrixPolicy Y(state.variables_); + MatrixPolicy Ynew(Y.size(), Y[0].size(), 0.0); + MatrixPolicy initial_forcing(Y.size(), Y[0].size(), 0.0); + MatrixPolicy forcing(Y.size(), Y[0].size(), 0.0); + MatrixPolicy temp(Y.size(), Y[0].size(), 0.0); + std::vector> K{}; + + parameters_.h_max_ = time_step; + parameters_.h_start_ = std::max(parameters_.h_min_, delta_min_); + + stats_.Reset(); + UpdateState(state); + + for (std::size_t i = 0; i < parameters_.stages_; ++i) + K.push_back(MatrixPolicy(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_)); + + if (std::abs(H) <= 10 * parameters_.round_off_) + H = delta_min_; + + bool reject_last_h = false; + bool reject_more_h = false; + + while ((present_time - time_step + parameters_.round_off_) <= 0 && (result.state_ == SolverState::Running)) + { + if (stats_.number_of_steps > parameters_.max_number_of_steps_) + { + result.state_ = SolverState::ConvergenceExceededMaxSteps; + break; + } + + if (((present_time + 0.1 * H) == present_time) || (H <= parameters_.round_off_)) + { + result.state_ = SolverState::StepSizeTooSmall; + break; + } + + // Limit H if necessary to avoid going beyond the specified chemistry time step + H = std::min(H, std::abs(time_step - present_time)); + + // compute the concentrations at the current time + CalculateForcing(state.rate_constants_, Y, initial_forcing); + + // compute the jacobian at the current time + CalculateJacobian(state.rate_constants_, Y, jacobian_); + + bool accepted = false; + // Repeat step calculation until current step accepted + while (!accepted) + { + bool is_singular{ false }; + // Form and factor the rosenbrock ode jacobian + LinearFactor(H, parameters_.gamma_[0], is_singular, Y, state.rate_constants_); + if (is_singular) + { + result.state_ = SolverState::RepeatedlySingularMatrix; + break; + } + + // Compute the stages + for (uint64_t stage = 0; stage < parameters_.stages_; ++stage) + { + double stage_combinations = ((stage + 1) - 1) * ((stage + 1) - 2) / 2; + if (stage == 0) + { + forcing = initial_forcing; + } + else + { + if (parameters_.new_function_evaluation_[stage]) + { + Ynew.AsVector().assign(Y.AsVector().begin(), Y.AsVector().end()); + for (uint64_t j = 0; j < stage; ++j) + { + auto a = parameters_.a_[stage_combinations + j]; + Ynew.ForEach([&](double& iYnew, double& iKj) { iYnew += a * iKj; }, K[j]); + } + CalculateForcing(state.rate_constants_, Ynew, forcing); + } + } + K[stage].AsVector().assign(forcing.AsVector().begin(), forcing.AsVector().end()); + for (uint64_t j = 0; j < stage; ++j) + { + auto HC = parameters_.c_[stage_combinations + j] / H; + K[stage].ForEach([&](double& iKstage, double& iKj) { iKstage += HC * iKj; }, K[j]); + } + temp.AsVector().assign(K[stage].AsVector().begin(), K[stage].AsVector().end()); + linear_solver_.template Solve(temp, K[stage]); + stats_.solves += 1; + } + + // Compute the new solution + Ynew.AsVector().assign(Y.AsVector().begin(), Y.AsVector().end()); + for (uint64_t stage = 0; stage < parameters_.stages_; ++stage) + Ynew.ForEach([&](double& iYnew, double& iKstage) { iYnew += parameters_.m_[stage] * iKstage; }, K[stage]); + + // Compute the error estimation + MatrixPolicy Yerror(Y.size(), Y[0].size(), 0); + for (uint64_t stage = 0; stage < parameters_.stages_; ++stage) + Yerror.ForEach([&](double& iYerror, double& iKstage) { iYerror += parameters_.e_[stage] * iKstage; }, K[stage]); + + auto error = NormalizedError(Y, Ynew, Yerror); + + // New step size is bounded by FacMin <= Hnew/H <= FacMax + double fac = std::min( + parameters_.factor_max_, + std::max( + parameters_.factor_min_, + parameters_.safety_factor_ / std::pow(error, 1 / parameters_.estimator_of_local_order_))); + double Hnew = H * fac; + + // Check the error magnitude and adjust step size + stats_.number_of_steps += 1; + stats_.total_steps += 1; + + if (std::isnan(error)) + { + Y.AsVector().assign(Ynew.AsVector().begin(), Ynew.AsVector().end()); + result.state_ = SolverState::NaNDetected; + break; + } + else if ((error < 1) || (H < parameters_.h_min_)) + { + 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_)); + if (reject_last_h) + { + // No step size increase after a rejected step + Hnew = std::min(Hnew, H); + } + reject_last_h = false; + reject_more_h = false; + H = Hnew; + accepted = true; + } + else + { + // Reject step + if (reject_more_h) + { + Hnew = H * parameters_.rejection_factor_decrease_; + } + reject_more_h = reject_last_h; + reject_last_h = true; + H = Hnew; + if (stats_.accepted >= 1) + { + stats_.rejected += 1; + } + } + } + } + + result.state_ = SolverState::Converged; + result.final_time_ = present_time; + result.stats_ = stats_; + result.result_ = std::move(Y); + return result; + } + + template class MatrixPolicy, template class SparseMatrixPolicy> + inline void RosenbrockSolver::CalculateForcing( + const MatrixPolicy& rate_constants, + const MatrixPolicy& number_densities, + MatrixPolicy& forcing) + { + std::fill(forcing.AsVector().begin(), forcing.AsVector().end(), 0.0); + process_set_.AddForcingTerms(rate_constants, number_densities, forcing); + stats_.function_calls += 1; + } + + template class MatrixPolicy, template class SparseMatrixPolicy> + inline void RosenbrockSolver::AlphaMinusJacobian( + SparseMatrixPolicy& jacobian, + const double& alpha) const + requires(!VectorizableSparse>) + { + for (auto& elem : jacobian.AsVector()) + elem = -elem; + for (std::size_t i_block = 0; i_block < jacobian.size(); ++i_block) + { + auto jacobian_vector = std::next(jacobian.AsVector().begin(), i_block * jacobian.FlatBlockSize()); + for (const auto& i_elem : jacobian_diagonal_elements_) + jacobian_vector[i_elem] += alpha; + } + } + + template class MatrixPolicy, template class SparseMatrixPolicy> + inline void RosenbrockSolver::AlphaMinusJacobian( + SparseMatrixPolicy& jacobian, + const double& alpha) const + requires(VectorizableSparse>) + { + const std::size_t n_cells = jacobian.GroupVectorSize(); + for (auto& elem : jacobian.AsVector()) + elem = -elem; + for (std::size_t i_group = 0; i_group < jacobian.NumberOfGroups(jacobian.size()); ++i_group) + { + auto jacobian_vector = std::next(jacobian.AsVector().begin(), i_group * jacobian.GroupSize(jacobian.FlatBlockSize())); + for (const auto& i_elem : jacobian_diagonal_elements_) + for (std::size_t i_cell = 0; i_cell < n_cells; ++i_cell) + jacobian_vector[i_elem + i_cell] += alpha; + } + } + + template class MatrixPolicy, template class SparseMatrixPolicy> + inline void RosenbrockSolver::CalculateJacobian( + const MatrixPolicy& rate_constants, + const MatrixPolicy& number_densities, + SparseMatrixPolicy& jacobian) + { + std::fill(jacobian.AsVector().begin(), jacobian.AsVector().end(), 0.0); + process_set_.AddJacobianTerms(rate_constants, number_densities, jacobian); + stats_.jacobian_updates += 1; + } + + template class MatrixPolicy, template class SparseMatrixPolicy> + inline void RosenbrockSolver::UpdateState(State& state) + { + Process::UpdateState(processes_, state); + } + + template class MatrixPolicy, template class SparseMatrixPolicy> + inline void RosenbrockSolver::LinearFactor( + double& H, + const double gamma, + bool& singular, + const MatrixPolicy& number_densities, + const MatrixPolicy& rate_constants) + { + // TODO: invesitage this function. The fortran equivalent appears to have a bug. + // From my understanding the fortran do loop would only ever do one iteration and is equivalent to what's below + SparseMatrixPolicy jacobian = jacobian_; + uint64_t n_consecutive = 0; + singular = true; + while (true) + { + double alpha = 1 / (H * gamma); + AlphaMinusJacobian(jacobian, alpha); + linear_solver_.Factor(jacobian); + singular = false; // TODO This should be evaluated in some way + stats_.decompositions += 1; + if (!singular) + break; + stats_.singular += 1; + if (++n_consecutive > 5) + break; + H /= 2; + } + } + + template class MatrixPolicy, template class SparseMatrixPolicy> + inline double RosenbrockSolver::NormalizedError( + const MatrixPolicy& Y, + const MatrixPolicy& Ynew, + const MatrixPolicy& errors) + { + // Solving Ordinary Differential Equations II, page 123 + // https://link-springer-com.cuucar.idm.oclc.org/book/10.1007/978-3-642-05221-7 + + auto _y = Y.AsVector(); + auto _ynew = Ynew.AsVector(); + auto _errors = errors.AsVector(); + + double error = 0; + + for (size_t i = 0; i < N_; ++i) + { + double ymax = std::max(std::abs(_y[i]), std::abs(_ynew[i])); + double scale = parameters_.absolute_tolerance_ + parameters_.relative_tolerance_ * ymax; + error += std::pow(_errors[i] / scale, 2); + } + + double error_min_ = 1.0e-10; + return std::max(std::sqrt(error / N_), error_min_); + } + +} // namespace micm \ No newline at end of file diff --git a/include/micm/solver/state.hpp b/include/micm/solver/state.hpp index e5e636077..abd56fcca 100644 --- a/include/micm/solver/state.hpp +++ b/include/micm/solver/state.hpp @@ -13,7 +13,7 @@ namespace micm { - + struct StateParameters { std::vector state_variable_names_{}; @@ -44,7 +44,7 @@ namespace micm /// @brief /// @param parameters State dimension information - State(const StateParameters parameters); + State(const StateParameters& parameters); /// @brief Set species' concentrations /// @param species_to_concentration @@ -68,157 +68,7 @@ namespace micm void SetCustomRateParameter(const std::string& label, double value); void SetCustomRateParameter(const std::string& label, const std::vector& values); }; + +} // namespace micm - template class MatrixPolicy> - inline State::State() - : conditions_(), - variable_map_(), - custom_rate_parameter_map_(), - variables_(), - custom_rate_parameters_(), - rate_constants_() - { - } - template class MatrixPolicy> - inline State::State( - const std::size_t state_size, - const std::size_t custom_parameters_size, - const std::size_t process_size) - : conditions_(1), - variable_map_(), - variable_names_(), - custom_rate_parameter_map_(), - variables_(1, state_size, 0.0), - custom_rate_parameters_(1, custom_parameters_size, 0.0), - rate_constants_(1, process_size, 0.0) - { - } - - template class MatrixPolicy> - inline State::State(const StateParameters parameters) - : conditions_(parameters.number_of_grid_cells_), - variable_map_(), - variable_names_(parameters.state_variable_names_), - custom_rate_parameter_map_(), - variables_(parameters.number_of_grid_cells_, parameters.state_variable_names_.size(), 0.0), - custom_rate_parameters_(parameters.number_of_grid_cells_, parameters.custom_rate_parameter_labels_.size(), 0.0), - rate_constants_(parameters.number_of_grid_cells_, parameters.number_of_rate_constants_, 0.0) - { - std::size_t index = 0; - for (auto& name : parameters.state_variable_names_) - variable_map_[name] = index++; - index = 0; - for (auto& label : parameters.custom_rate_parameter_labels_) - custom_rate_parameter_map_[label] = index++; - } - - template class MatrixPolicy> - inline void State::SetConcentrations( - const System& system, - const std::unordered_map>& species_to_concentration) - { - int num_set_grid_cells = 0; - unsigned num_species = system.gas_phase_.species_.size(); - - std::vector num_concentrations_per_species; - num_concentrations_per_species.reserve(num_species); - - // Iterate map to store the number of concentration values corresponding to the number of set of grid cells - for (auto& species : system.gas_phase_.species_) - { - auto species_ptr = species_to_concentration.find(species.name_); - if (species_ptr == species_to_concentration.end()) - { - throw std::runtime_error("Concentration value(s) for '" + species.name_ + "' must be given."); - } - num_concentrations_per_species.push_back(species_ptr->second.size()); - } - - // Check if number of concentraiton inputs are the same for all species - if (!std::all_of( - num_concentrations_per_species.begin(), - num_concentrations_per_species.end(), - [&](int& i) { return i == num_concentrations_per_species.front(); })) - { - throw std::runtime_error("Concentration value must be given to all sets of grid cells."); - } - - num_set_grid_cells = num_concentrations_per_species[0]; - - // Find species and iterate through the keys to store concentrations for each set of grid cells - // 'concentrations' represents an N-D array in contiguous memory (N = num_set_grid_cells) - std::vector concentrations; - concentrations.resize(num_species * num_set_grid_cells); - - for (int i = 0; i < num_species; i++) - { - auto species_ptr = species_to_concentration.find(system.gas_phase_.species_[i].name_); - - for (int j = 0; j < num_set_grid_cells; j++) - { - concentrations[i + num_species * j] = species_ptr->second[j]; - } - } - - // Extract sub vector to assign to the corresponding set of grid cells. - // TODO: jiwon 7/12 - I think we want to reduce copy operations here - std::vector sub_concentrations; - sub_concentrations.reserve(num_species); - - for (int i = 0; i < num_set_grid_cells; i++) - { - sub_concentrations = { concentrations.begin() + (num_species * i), - concentrations.begin() + (num_species * i) + num_species }; - variables_[i] = sub_concentrations; - } - } - - template class MatrixPolicy> - void State::SetConcentration(const Species& species, double concentration) - { - if (variables_.size() != 1) - throw std::invalid_argument("Incorrect number of concentration values passed to multi-gridcell State"); - variables_[0][variable_map_[species.name_]] = concentration; - } - - template class MatrixPolicy> - void State::SetConcentration(const Species& species, const std::vector& concentration) - { - if (variables_.size() != concentration.size()) - throw std::invalid_argument("Incorrect number of concentration values passed to multi-gridcell State"); - std::size_t i_species = variable_map_[species.name_]; - for (std::size_t i = 0; i < variables_.size(); ++i) - variables_[i][i_species] = concentration[i]; - } - - template class MatrixPolicy> - inline void State::SetCustomRateParameters( - const std::unordered_map>& parameters) - { - for (auto& pair : parameters) - SetCustomRateParameter(pair.first, pair.second); - } - - template class MatrixPolicy> - inline void State::SetCustomRateParameter(const std::string& label, double value) - { - auto param = custom_rate_parameter_map_.find(label); - if (param == custom_rate_parameter_map_.end()) - throw std::invalid_argument("Unkonwn rate constant parameter '" + label + "'"); - if (custom_rate_parameters_.size() != 1) - throw std::invalid_argument("Incorrect number of custom rate parameter values passed to multi-gridcell State"); - custom_rate_parameters_[0][param->second] = value; - } - - template class MatrixPolicy> - inline void State::SetCustomRateParameter(const std::string& label, const std::vector& values) - { - auto param = custom_rate_parameter_map_.find(label); - if (param == custom_rate_parameter_map_.end()) - throw std::invalid_argument("Unkonwn rate constant parameter '" + label + "'"); - if (custom_rate_parameters_.size() != values.size()) - throw std::invalid_argument("Incorrect number of custom rate parameter values passed to multi-gridcell State"); - for (std::size_t i = 0; i < custom_rate_parameters_.size(); ++i) - custom_rate_parameters_[i][param->second] = values[i]; - } -} // namespace micm \ No newline at end of file +#include "state.inl" \ No newline at end of file diff --git a/include/micm/solver/state.inl b/include/micm/solver/state.inl new file mode 100644 index 000000000..4fdfb8c0e --- /dev/null +++ b/include/micm/solver/state.inl @@ -0,0 +1,157 @@ +namespace micm +{ + + template class MatrixPolicy> + inline State::State() + : conditions_(), + variable_map_(), + custom_rate_parameter_map_(), + variable_names_(), + variables_(), + custom_rate_parameters_(), + rate_constants_() + { + } + + template class MatrixPolicy> + inline State::State( + const std::size_t state_size, + const std::size_t custom_parameters_size, + const std::size_t process_size) + : conditions_(1), + variable_map_(), + custom_rate_parameter_map_(), + variable_names_(), + variables_(1, state_size, 0.0), + custom_rate_parameters_(1, custom_parameters_size, 0.0), + rate_constants_(1, process_size, 0.0) + { + } + + template class MatrixPolicy> + inline State::State(const StateParameters& parameters) + : conditions_(parameters.number_of_grid_cells_), + variable_map_(), + custom_rate_parameter_map_(), + variable_names_(parameters.state_variable_names_), + variables_(parameters.number_of_grid_cells_, parameters.state_variable_names_.size(), 0.0), + custom_rate_parameters_(parameters.number_of_grid_cells_, parameters.custom_rate_parameter_labels_.size(), 0.0), + rate_constants_(parameters.number_of_grid_cells_, parameters.number_of_rate_constants_, 0.0) + { + std::size_t index = 0; + for (auto& name : parameters.state_variable_names_) + variable_map_[name] = index++; + index = 0; + for (auto& label : parameters.custom_rate_parameter_labels_) + custom_rate_parameter_map_[label] = index++; + } + + template class MatrixPolicy> + inline void State::SetConcentrations( + const System& system, + const std::unordered_map>& species_to_concentration) + { + int num_set_grid_cells = 0; + unsigned num_species = system.gas_phase_.species_.size(); + + std::vector num_concentrations_per_species; + num_concentrations_per_species.reserve(num_species); + + // Iterate map to store the number of concentration values corresponding to the number of set of grid cells + for (auto& species : system.gas_phase_.species_) + { + auto species_ptr = species_to_concentration.find(species.name_); + if (species_ptr == species_to_concentration.end()) + { + throw std::runtime_error("Concentration value(s) for '" + species.name_ + "' must be given."); + } + num_concentrations_per_species.push_back(species_ptr->second.size()); + } + + // Check if number of concentraiton inputs are the same for all species + if (!std::all_of( + num_concentrations_per_species.begin(), + num_concentrations_per_species.end(), + [&](int& i) { return i == num_concentrations_per_species.front(); })) + { + throw std::runtime_error("Concentration value must be given to all sets of grid cells."); + } + + num_set_grid_cells = num_concentrations_per_species[0]; + + // Find species and iterate through the keys to store concentrations for each set of grid cells + // 'concentrations' represents an N-D array in contiguous memory (N = num_set_grid_cells) + std::vector concentrations; + concentrations.resize(num_species * num_set_grid_cells); + + for (int i = 0; i < num_species; i++) + { + auto species_ptr = species_to_concentration.find(system.gas_phase_.species_[i].name_); + + for (int j = 0; j < num_set_grid_cells; j++) + { + concentrations[i + num_species * j] = species_ptr->second[j]; + } + } + + // Extract sub vector to assign to the corresponding set of grid cells. + std::vector sub_concentrations; + sub_concentrations.reserve(num_species); + + for (int i = 0; i < num_set_grid_cells; i++) + { + sub_concentrations = { concentrations.begin() + (num_species * i), + concentrations.begin() + (num_species * i) + num_species }; + variables_[i] = sub_concentrations; + } + } + + template class MatrixPolicy> + inline void State::SetConcentration(const Species& species, double concentration) + { + if (variables_.size() != 1) + throw std::invalid_argument("Incorrect number of concentration values passed to multi-gridcell State"); + variables_[0][variable_map_[species.name_]] = concentration; + } + + template class MatrixPolicy> + inline void State::SetConcentration(const Species& species, const std::vector& concentration) + { + if (variables_.size() != concentration.size()) + throw std::invalid_argument("Incorrect number of concentration values passed to multi-gridcell State"); + std::size_t i_species = variable_map_[species.name_]; + for (std::size_t i = 0; i < variables_.size(); ++i) + variables_[i][i_species] = concentration[i]; + } + + template class MatrixPolicy> + inline void State::SetCustomRateParameters(const std::unordered_map>& parameters) + { + for (auto& pair : parameters) + SetCustomRateParameter(pair.first, pair.second); + } + + template class MatrixPolicy> + inline void State::SetCustomRateParameter(const std::string& label, double value) + { + auto param = custom_rate_parameter_map_.find(label); + if (param == custom_rate_parameter_map_.end()) + throw std::invalid_argument("Unkonwn rate constant parameter '" + label + "'"); + if (custom_rate_parameters_.size() != 1) + throw std::invalid_argument("Incorrect number of custom rate parameter values passed to multi-gridcell State"); + custom_rate_parameters_[0][param->second] = value; + } + + template class MatrixPolicy> + inline void State::SetCustomRateParameter(const std::string& label, const std::vector& values) + { + auto param = custom_rate_parameter_map_.find(label); + if (param == custom_rate_parameter_map_.end()) + throw std::invalid_argument("Unkonwn rate constant parameter '" + label + "'"); + if (custom_rate_parameters_.size() != values.size()) + throw std::invalid_argument("Incorrect number of custom rate parameter values passed to multi-gridcell State"); + for (std::size_t i = 0; i < custom_rate_parameters_.size(); ++i) + custom_rate_parameters_[i][param->second] = values[i]; + } + +} // namespace micm