diff --git a/include/micm/process/cuda_process_set.hpp b/include/micm/process/cuda_process_set.hpp index f678ac30d..ef72f1c01 100644 --- a/include/micm/process/cuda_process_set.hpp +++ b/include/micm/process/cuda_process_set.hpp @@ -102,11 +102,10 @@ namespace micm } template - requires(CudaMatrix&& VectorizableDense) inline void CudaProcessSet:: - AddForcingTerms( - const MatrixPolicy& rate_constants, - const MatrixPolicy& state_variables, - MatrixPolicy& forcing) const + requires(CudaMatrix&& VectorizableDense) inline void CudaProcessSet::AddForcingTerms( + const MatrixPolicy& rate_constants, + const MatrixPolicy& state_variables, + MatrixPolicy& forcing) const { auto forcing_param = forcing.AsDeviceParam(); // we need to update forcing so it can't be constant and must be an lvalue micm::cuda::AddForcingTermsKernelDriver( diff --git a/include/micm/process/jit_process_set.hpp b/include/micm/process/jit_process_set.hpp index aaeabbd13..d5ab35145 100644 --- a/include/micm/process/jit_process_set.hpp +++ b/include/micm/process/jit_process_set.hpp @@ -34,9 +34,7 @@ namespace micm /// @brief Create a JITed process set calculator for a given set of processes /// @param processes Processes to create calculator for /// @param variable_map A mapping of species names to concentration index - JitProcessSet( - const std::vector &processes, - const std::map &variable_map); + JitProcessSet(const std::vector &processes, const std::map &variable_map); ~JitProcessSet(); @@ -50,10 +48,8 @@ namespace micm /// @param state_variables Current state variable values (grid cell, state variable) /// @param forcing Forcing terms for each state variable (grid cell, state variable) template - void AddForcingTerms( - const MatrixPolicy &rate_constants, - const MatrixPolicy &state_variables, - MatrixPolicy &forcing) const; + void AddForcingTerms(const MatrixPolicy &rate_constants, const MatrixPolicy &state_variables, MatrixPolicy &forcing) + const; /// @brief Subtracts Jacobian terms for the set of processes for the current conditions /// @param rate_constants Current values for the process rate constants (grid cell, process) diff --git a/include/micm/solver/backward_euler.hpp b/include/micm/solver/backward_euler.hpp index 8288a51bd..b5d1fb404 100644 --- a/include/micm/solver/backward_euler.hpp +++ b/include/micm/solver/backward_euler.hpp @@ -29,7 +29,7 @@ namespace micm { /// @brief An implementation of the fully implicit backward euler method - template + template class BackwardEuler { BackwardEulerSolverParameters parameters_; @@ -41,12 +41,12 @@ namespace micm public: /// @brief Default constructor BackwardEuler( - BackwardEulerSolverParameters parameters, - LinearSolverPolicy linear_solver, - ProcessSetPolicy process_set, - std::vector jacobian_diagonal_elements, - std::vector& processes - ) : parameters_(parameters), + BackwardEulerSolverParameters parameters, + LinearSolverPolicy linear_solver, + ProcessSetPolicy process_set, + std::vector jacobian_diagonal_elements, + std::vector& processes) + : parameters_(parameters), linear_solver_(linear_solver), process_set_(process_set), jacobian_diagonal_elements_(jacobian_diagonal_elements), @@ -60,9 +60,7 @@ namespace micm /// @param time_step Time [s] to advance the state by /// @param state The state to advance /// @return Nothing, but the state is updated - void Solve( - double time_step, - auto& state); + void Solve(double time_step, auto& state); }; } // namespace micm diff --git a/include/micm/solver/backward_euler.inl b/include/micm/solver/backward_euler.inl index 7bae79d9f..f51d4fcf0 100644 --- a/include/micm/solver/backward_euler.inl +++ b/include/micm/solver/backward_euler.inl @@ -44,10 +44,8 @@ inline std::error_code make_error_code(MicmBackwardEulerErrc e) namespace micm { - template - inline void BackwardEuler::Solve( - double time_step, - auto& state) + template + inline void BackwardEuler::Solve(double time_step, auto& state) { // A fully implicit euler implementation is given by the following equation: // y_{n+1} = y_n + H * f(t_{n+1}, y_{n+1}) @@ -60,7 +58,6 @@ namespace micm // if that fails, try H = H/10 once // if that fails, accept the current H but do not update the Yn vector - double tolerance = parameters_.absolute_tolerance_[0]; double small = parameters_.small; std::size_t max_iter = parameters_.max_number_of_steps_; diff --git a/include/micm/solver/backward_euler_solver_parameters.hpp b/include/micm/solver/backward_euler_solver_parameters.hpp index 5aa211c0a..a1fc00aad 100644 --- a/include/micm/solver/backward_euler_solver_parameters.hpp +++ b/include/micm/solver/backward_euler_solver_parameters.hpp @@ -4,10 +4,10 @@ */ #pragma once +#include #include #include #include -#include namespace micm { @@ -16,7 +16,7 @@ namespace micm struct BackwardEulerSolverParameters { std::vector absolute_tolerance_; - double small { 1.0e-40 }; + double small{ 1.0e-40 }; size_t max_number_of_steps_{ 11 }; std::array time_step_reductions{ 0.5, 0.5, 0.5, 0.5, 0.1 }; }; diff --git a/include/micm/solver/cuda_linear_solver.hpp b/include/micm/solver/cuda_linear_solver.hpp index bba93ab83..0b36e0c91 100644 --- a/include/micm/solver/cuda_linear_solver.hpp +++ b/include/micm/solver/cuda_linear_solver.hpp @@ -70,8 +70,8 @@ namespace micm template requires( - CudaMatrix&& CudaMatrix&& VectorizableDense&& - VectorizableSparse) void Solve(const MatrixPolicy& b, MatrixPolicy& x, const SparseMatrixPolicy& L, const SparseMatrixPolicy& U) + CudaMatrix&& CudaMatrix&& VectorizableDense&& VectorizableSparse< + SparseMatrixPolicy>) void Solve(const MatrixPolicy& b, MatrixPolicy& x, const SparseMatrixPolicy& L, const SparseMatrixPolicy& U) const { auto x_param = x.AsDeviceParam(); // we need to update x so it can't be constant and must be an lvalue diff --git a/include/micm/solver/cuda_lu_decomposition.hpp b/include/micm/solver/cuda_lu_decomposition.hpp index 38e659f92..83010d91a 100644 --- a/include/micm/solver/cuda_lu_decomposition.hpp +++ b/include/micm/solver/cuda_lu_decomposition.hpp @@ -89,8 +89,10 @@ namespace micm }; template - requires(CudaMatrix&& VectorizableSparse) void CudaLuDecomposition:: - Decompose(const SparseMatrixPolicy& A, SparseMatrixPolicy& L, SparseMatrixPolicy& U) const + requires(CudaMatrix&& VectorizableSparse) void CudaLuDecomposition::Decompose( + const SparseMatrixPolicy& A, + SparseMatrixPolicy& L, + SparseMatrixPolicy& U) const { auto L_param = L.AsDeviceParam(); // we need to update lower matrix so it can't be constant and must be an lvalue auto U_param = U.AsDeviceParam(); // we need to update upper matrix so it can't be constant and must be an lvalue diff --git a/include/micm/solver/jit_linear_solver.hpp b/include/micm/solver/jit_linear_solver.hpp index 0c5d55a56..5a660dfe6 100644 --- a/include/micm/solver/jit_linear_solver.hpp +++ b/include/micm/solver/jit_linear_solver.hpp @@ -37,9 +37,7 @@ namespace micm /// @param compiler JIT compiler /// @param matrix Block-diagonal sparse matrix to create solver for /// @param initial_value Initial value for decomposed triangular matrix elements - JitLinearSolver( - const SparseMatrixPolicy& matrix, - double initial_value); + JitLinearSolver(const SparseMatrixPolicy& matrix, double initial_value); ~JitLinearSolver(); @@ -54,18 +52,11 @@ namespace micm /// @brief Decompose the matrix into upper and lower triangular matrices and general JIT functions /// @param matrix Matrix that will be factored into lower and upper triangular matrices - void Factor( - SparseMatrixPolicy& matrix, - SparseMatrixPolicy& lower_matrix, - SparseMatrixPolicy& upper_matrix); + void Factor(SparseMatrixPolicy& matrix, SparseMatrixPolicy& lower_matrix, SparseMatrixPolicy& upper_matrix); /// @brief Solve for x in Ax = b template - void Solve( - const MatrixPolicy& b, - MatrixPolicy& x, - SparseMatrixPolicy& lower_matrix, - SparseMatrixPolicy& upper_matrix); + void Solve(const MatrixPolicy& b, MatrixPolicy& x, SparseMatrixPolicy& lower_matrix, SparseMatrixPolicy& upper_matrix); private: /// @brief Generates the JIT-ed Solve function diff --git a/include/micm/solver/jit_linear_solver.inl b/include/micm/solver/jit_linear_solver.inl index ff40277eb..bb066902f 100644 --- a/include/micm/solver/jit_linear_solver.inl +++ b/include/micm/solver/jit_linear_solver.inl @@ -32,8 +32,7 @@ namespace micm : LinearSolver( matrix, initial_value, - [&](const SparseMatrixPolicy &m) -> LuDecompositionPolicy - { return LuDecompositionPolicy(m); }) + [&](const SparseMatrixPolicy &m) -> LuDecompositionPolicy { return LuDecompositionPolicy(m); }) { solve_function_ = NULL; if (matrix.NumberOfBlocks() != L || matrix.GroupVectorSize() != L) diff --git a/include/micm/solver/jit_lu_decomposition.hpp b/include/micm/solver/jit_lu_decomposition.hpp index 08a01e091..0f24ddd2d 100644 --- a/include/micm/solver/jit_lu_decomposition.hpp +++ b/include/micm/solver/jit_lu_decomposition.hpp @@ -42,11 +42,8 @@ namespace micm /// @param upper The upper triangular matrix created by decomposition /// @param is_singular Flag that will be set to true if A is singular; false otherwise template - void Decompose( - const SparseMatrixPolicy &A, - SparseMatrixPolicy &lower, - SparseMatrixPolicy &upper, - bool &is_singular) const; + void Decompose(const SparseMatrixPolicy &A, SparseMatrixPolicy &lower, SparseMatrixPolicy &upper, bool &is_singular) + const; /// @brief Create sparse L and U matrices for a given A matrix /// @param A Sparse matrix that will be decomposed diff --git a/include/micm/solver/jit_lu_decomposition.inl b/include/micm/solver/jit_lu_decomposition.inl index a1004c864..c3aafe93b 100644 --- a/include/micm/solver/jit_lu_decomposition.inl +++ b/include/micm/solver/jit_lu_decomposition.inl @@ -25,8 +25,7 @@ namespace micm } template - inline JitLuDecomposition::JitLuDecomposition( - const SparseMatrix> &matrix) + inline JitLuDecomposition::JitLuDecomposition(const SparseMatrix> &matrix) : LuDecomposition(LuDecomposition::Create>>(matrix)) { decompose_function_ = NULL; @@ -198,10 +197,8 @@ namespace micm template template - void JitLuDecomposition::Decompose( - const SparseMatrixPolicy &A, - SparseMatrixPolicy &lower, - SparseMatrixPolicy &upper) const + void JitLuDecomposition::Decompose(const SparseMatrixPolicy &A, SparseMatrixPolicy &lower, SparseMatrixPolicy &upper) + const { decompose_function_(A.AsVector().data(), lower.AsVector().data(), upper.AsVector().data()); } diff --git a/include/micm/solver/linear_solver.inl b/include/micm/solver/linear_solver.inl index d252b9aa6..ecd304b5f 100644 --- a/include/micm/solver/linear_solver.inl +++ b/include/micm/solver/linear_solver.inl @@ -136,7 +136,8 @@ namespace micm template requires( !VectorizableDense || - !VectorizableSparse) inline void LinearSolver::Solve( + !VectorizableSparse) inline void LinearSolver:: + Solve( const MatrixPolicy& b, MatrixPolicy& x, const SparseMatrixPolicy& lower_matrix, @@ -191,8 +192,8 @@ namespace micm template template - requires(VectorizableDense&& VectorizableSparse< - SparseMatrixPolicy>) inline void LinearSolver:: + requires(VectorizableDense&& + VectorizableSparse) inline void LinearSolver:: Solve( const MatrixPolicy& b, MatrixPolicy& x, diff --git a/include/micm/solver/lu_decomposition.hpp b/include/micm/solver/lu_decomposition.hpp index 3173f1053..538818a79 100644 --- a/include/micm/solver/lu_decomposition.hpp +++ b/include/micm/solver/lu_decomposition.hpp @@ -90,7 +90,9 @@ namespace micm /// @param A Sparse matrix that will be decomposed /// @return L and U Sparse matrices template - static std::pair GetLUMatrices(const SparseMatrixPolicy& A, typename SparseMatrixPolicy::value_type initial_value); + static std::pair GetLUMatrices( + const SparseMatrixPolicy& A, + typename SparseMatrixPolicy::value_type initial_value); /// @brief Perform an LU decomposition on a given A matrix /// @param A Sparse matrix to decompose diff --git a/include/micm/solver/lu_decomposition.inl b/include/micm/solver/lu_decomposition.inl index d9803d42b..dc378db4f 100644 --- a/include/micm/solver/lu_decomposition.inl +++ b/include/micm/solver/lu_decomposition.inl @@ -162,8 +162,7 @@ namespace micm } template - inline void LuDecomposition::Decompose(const SparseMatrixPolicy& A, SparseMatrixPolicy& L, SparseMatrixPolicy& U) - const + inline void LuDecomposition::Decompose(const SparseMatrixPolicy& A, SparseMatrixPolicy& L, SparseMatrixPolicy& U) const { bool is_singular = false; Decompose(A, L, U, is_singular); diff --git a/include/micm/solver/rosenbrock.inl b/include/micm/solver/rosenbrock.inl index 9a48b3c07..3c9ba5676 100644 --- a/include/micm/solver/rosenbrock.inl +++ b/include/micm/solver/rosenbrock.inl @@ -75,12 +75,7 @@ namespace micm } } - template< - template - class MatrixPolicy, - class SparseMatrixPolicy, - class LinearSolverPolicy, - class ProcessSetPolicy> + template class MatrixPolicy, class SparseMatrixPolicy, class LinearSolverPolicy, class ProcessSetPolicy> inline RosenbrockSolver::RosenbrockSolver() : processes_(), parameters_(RosenbrockSolverParameters::ThreeStageRosenbrockParameters()), @@ -90,12 +85,7 @@ namespace micm { } - template< - template - class MatrixPolicy, - class SparseMatrixPolicy, - class LinearSolverPolicy, - class ProcessSetPolicy> + template class MatrixPolicy, class SparseMatrixPolicy, class LinearSolverPolicy, class ProcessSetPolicy> inline RosenbrockSolver::RosenbrockSolver( const System& system, const std::vector& processes, @@ -114,12 +104,7 @@ namespace micm { } - template< - template - class MatrixPolicy, - class SparseMatrixPolicy, - class LinearSolverPolicy, - class ProcessSetPolicy> + template class MatrixPolicy, class SparseMatrixPolicy, class LinearSolverPolicy, class ProcessSetPolicy> inline RosenbrockSolver::RosenbrockSolver( const System& system, const std::vector& processes, @@ -231,24 +216,14 @@ namespace micm linear_solver_ = std::move(create_linear_solver(jacobian, 1.0e-30)); } - template< - template - class MatrixPolicy, - class SparseMatrixPolicy, - class LinearSolverPolicy, - class ProcessSetPolicy> + template class MatrixPolicy, class SparseMatrixPolicy, class LinearSolverPolicy, class ProcessSetPolicy> inline State, SparseMatrixPolicy> RosenbrockSolver::GetState() const { return State, SparseMatrixPolicy>{ state_parameters_ }; } - template< - template - class MatrixPolicy, - class SparseMatrixPolicy, - class LinearSolverPolicy, - class ProcessSetPolicy> + template class MatrixPolicy, class SparseMatrixPolicy, class LinearSolverPolicy, class ProcessSetPolicy> inline typename RosenbrockSolver::SolverResult RosenbrockSolver::Solve( double time_step, @@ -447,12 +422,7 @@ namespace micm return result; } - template< - template - class MatrixPolicy, - class SparseMatrixPolicy, - class LinearSolverPolicy, - class ProcessSetPolicy> + template class MatrixPolicy, class SparseMatrixPolicy, class LinearSolverPolicy, class ProcessSetPolicy> inline void RosenbrockSolver::CalculateForcing( const MatrixPolicy& rate_constants, const MatrixPolicy& number_densities, @@ -463,12 +433,7 @@ namespace micm process_set_.template AddForcingTerms>(rate_constants, number_densities, forcing); } - template< - template - class MatrixPolicy, - class SparseMatrixPolicy, - class LinearSolverPolicy, - class ProcessSetPolicy> + template class MatrixPolicy, class SparseMatrixPolicy, class LinearSolverPolicy, class ProcessSetPolicy> inline void RosenbrockSolver::AlphaMinusJacobian( SparseMatrixPolicy& jacobian, const double& alpha) const requires(!VectorizableSparse) @@ -483,12 +448,7 @@ namespace micm } } - template< - template - class MatrixPolicy, - class SparseMatrixPolicy, - class LinearSolverPolicy, - class ProcessSetPolicy> + template class MatrixPolicy, class SparseMatrixPolicy, class LinearSolverPolicy, class ProcessSetPolicy> inline void RosenbrockSolver::AlphaMinusJacobian( SparseMatrixPolicy& jacobian, const double& alpha) const requires(VectorizableSparse) @@ -505,12 +465,7 @@ namespace micm } } - template< - template - class MatrixPolicy, - class SparseMatrixPolicy, - class LinearSolverPolicy, - class ProcessSetPolicy> + template class MatrixPolicy, class SparseMatrixPolicy, class LinearSolverPolicy, class ProcessSetPolicy> inline void RosenbrockSolver::CalculateNegativeJacobian( const MatrixPolicy& rate_constants, @@ -523,24 +478,14 @@ namespace micm process_set_.SubtractJacobianTerms(rate_constants, number_densities, jacobian); } - template< - template - class MatrixPolicy, - class SparseMatrixPolicy, - class LinearSolverPolicy, - class ProcessSetPolicy> + template class MatrixPolicy, class SparseMatrixPolicy, class LinearSolverPolicy, class ProcessSetPolicy> inline void RosenbrockSolver::UpdateState( State, SparseMatrixPolicy>& state) { Process::UpdateState(processes_, state); } - template< - template - class MatrixPolicy, - class SparseMatrixPolicy, - class LinearSolverPolicy, - class ProcessSetPolicy> + template class MatrixPolicy, class SparseMatrixPolicy, class LinearSolverPolicy, class ProcessSetPolicy> inline void RosenbrockSolver::LinearFactor( double& H, const double gamma, @@ -578,12 +523,7 @@ namespace micm } } - template< - template - class MatrixPolicy, - class SparseMatrixPolicy, - class LinearSolverPolicy, - class ProcessSetPolicy> + template class MatrixPolicy, class SparseMatrixPolicy, class LinearSolverPolicy, class ProcessSetPolicy> inline double RosenbrockSolver::NormalizedError( const MatrixPolicy& Y, const MatrixPolicy& Ynew, @@ -617,12 +557,7 @@ namespace micm return std::max(std::sqrt(error / N), error_min); } - template< - template - class MatrixPolicy, - class SparseMatrixPolicy, - class LinearSolverPolicy, - class ProcessSetPolicy> + template class MatrixPolicy, class SparseMatrixPolicy, class LinearSolverPolicy, class ProcessSetPolicy> inline double RosenbrockSolver::NormalizedError( const MatrixPolicy& Y, const MatrixPolicy& Ynew, diff --git a/include/micm/solver/rosenbrock_solver_parameters.hpp b/include/micm/solver/rosenbrock_solver_parameters.hpp index 46fc79c12..73d4d28cd 100644 --- a/include/micm/solver/rosenbrock_solver_parameters.hpp +++ b/include/micm/solver/rosenbrock_solver_parameters.hpp @@ -59,10 +59,10 @@ namespace micm std::vector absolute_tolerance_; double relative_tolerance_{ 1e-4 }; - std::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 - bool check_singularity_{ false }; // Check for singular A matrix in linear solve of A x = b - bool ignore_unused_species_{ false }; // Allow unused species to be included in state and solve + std::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 + bool check_singularity_{ false }; // Check for singular A matrix in linear solve of A x = b + bool ignore_unused_species_{ false }; // Allow unused species to be included in state and solve // Print RosenbrockSolverParameters to console void Print() const; diff --git a/include/micm/solver/solver.hpp b/include/micm/solver/solver.hpp index 33dea939c..382c2c330 100644 --- a/include/micm/solver/solver.hpp +++ b/include/micm/solver/solver.hpp @@ -4,7 +4,6 @@ */ #pragma once - namespace micm { @@ -19,7 +18,6 @@ namespace micm SolverPolicy solver_; public: - Solver( SolverPolicy solver, StateParameters state_parameters, diff --git a/include/micm/solver/solver_builder.hpp b/include/micm/solver/solver_builder.hpp index 5b91abb69..ab9d81152 100644 --- a/include/micm/solver/solver_builder.hpp +++ b/include/micm/solver/solver_builder.hpp @@ -18,8 +18,8 @@ #include #include -#include #include +#include namespace micm { @@ -69,7 +69,8 @@ namespace micm protected: /// @brief /// @return - virtual Solver, ProcessSet>, State> BuildBackwardEulerSolver() = 0; + virtual Solver, ProcessSet>, State> + BuildBackwardEulerSolver() = 0; virtual ~SolverBuilder() = default; @@ -101,14 +102,16 @@ namespace micm class CpuSolverBuilder : public SolverBuilder { public: - Solver, ProcessSet>, State> BuildBackwardEulerSolver() override; + Solver, ProcessSet>, State> + BuildBackwardEulerSolver() override; }; template class GpuSolverBuilder : public SolverBuilder { public: - Solver, ProcessSet>, State> BuildBackwardEulerSolver() override; + Solver, ProcessSet>, State> + BuildBackwardEulerSolver() override; }; } // namespace micm diff --git a/include/micm/solver/solver_builder.inl b/include/micm/solver/solver_builder.inl index 504d8d9d1..0116b2967 100644 --- a/include/micm/solver/solver_builder.inl +++ b/include/micm/solver/solver_builder.inl @@ -229,9 +229,7 @@ namespace micm } template - inline Solver< - BackwardEuler, ProcessSet>, - State> + inline Solver, ProcessSet>, State> CpuSolverBuilder::BuildBackwardEulerSolver() { using ProcessSetPolicy = ProcessSet; diff --git a/include/micm/solver/state.inl b/include/micm/solver/state.inl index 186b7d0f8..447bb0cbb 100644 --- a/include/micm/solver/state.inl +++ b/include/micm/solver/state.inl @@ -58,7 +58,7 @@ inline std::error_code make_error_code(MicmStateErrc e) namespace micm { - template + template inline State::State() : conditions_(), variables_(), diff --git a/include/micm/util/cuda_dense_matrix.hpp b/include/micm/util/cuda_dense_matrix.hpp index 41a5c517a..d36aca792 100644 --- a/include/micm/util/cuda_dense_matrix.hpp +++ b/include/micm/util/cuda_dense_matrix.hpp @@ -62,7 +62,7 @@ namespace micm template class CudaDenseMatrix : public VectorMatrix { - public: + public: // Diagonal markowitz reordering requires an int argument, make sure one is always accessible using IntMatrix = CudaDenseMatrix; using value_type = T; diff --git a/include/micm/util/cuda_sparse_matrix.hpp b/include/micm/util/cuda_sparse_matrix.hpp index aa0a6bfd4..bacc4352b 100644 --- a/include/micm/util/cuda_sparse_matrix.hpp +++ b/include/micm/util/cuda_sparse_matrix.hpp @@ -20,7 +20,7 @@ namespace micm template class CudaSparseMatrix : public SparseMatrix { - public: + public: // Diagonal markowitz reordering requires an int argument, make sure one is always accessible using IntMatrix = CudaSparseMatrix; using value_type = T; diff --git a/include/micm/util/matrix.hpp b/include/micm/util/matrix.hpp index 0b2b4c2ee..7ba008fca 100644 --- a/include/micm/util/matrix.hpp +++ b/include/micm/util/matrix.hpp @@ -28,12 +28,12 @@ namespace micm template class Matrix { - public: + public: // Diagonal markowitz reordering requires an int argument, make sure one is always accessible using IntMatrix = Matrix; using value_type = T; - private: + private: std::vector data_; std::size_t x_dim_; std::size_t y_dim_; diff --git a/include/micm/util/vector_matrix.hpp b/include/micm/util/vector_matrix.hpp index 12a368024..2a3bafd8c 100644 --- a/include/micm/util/vector_matrix.hpp +++ b/include/micm/util/vector_matrix.hpp @@ -31,12 +31,12 @@ namespace micm template class VectorMatrix { - public: + public: // Diagonal markowitz reordering requires an int argument, make sure one is always accessible using IntMatrix = VectorMatrix; using value_type = T; - private: + private: protected: std::vector data_; std::size_t x_dim_; // number of rows diff --git a/test/integration/analytical_policy.hpp b/test/integration/analytical_policy.hpp index a44a8202b..abb0b559d 100644 --- a/test/integration/analytical_policy.hpp +++ b/test/integration/analytical_policy.hpp @@ -110,7 +110,8 @@ void test_analytical_troe( .SetPhase(gas_phase); auto processes = std::vector{ r1, r2 }; - OdeSolverPolicy solver = OdeSolverPolicy(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), processes, parameters); + OdeSolverPolicy solver = + OdeSolverPolicy(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), processes, parameters); auto be_state = solver.GetState(); @@ -256,7 +257,9 @@ void test_analytical_stiff_troe( .SetPhase(gas_phase); OdeSolverPolicy solver = OdeSolverPolicy( - micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ r1, r2, r3, r4, r5 }, parameters); + micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), + std::vector{ r1, r2, r3, r4, r5 }, + parameters); double temperature = 272.5; double pressure = 101253.3; @@ -365,8 +368,8 @@ void test_analytical_photolysis( .SetRateConstant(micm::UserDefinedRateConstant({ .label_ = "photoB" })) .SetPhase(gas_phase); - OdeSolverPolicy solver = - OdeSolverPolicy(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ r1, r2 }, parameters); + OdeSolverPolicy solver = OdeSolverPolicy( + micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ r1, r2 }, parameters); double temperature = 272.5; double pressure = 101253.3; @@ -493,7 +496,9 @@ void test_analytical_stiff_photolysis( .SetPhase(gas_phase); OdeSolverPolicy solver = OdeSolverPolicy( - micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ r1, r2, r3, r4, r5 }, parameters); + micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), + std::vector{ r1, r2, r3, r4, r5 }, + parameters); double temperature = 272.5; double pressure = 101253.3; @@ -609,8 +614,8 @@ void test_analytical_ternary_chemical_activation( .N_ = 0.8 })) .SetPhase(gas_phase); - OdeSolverPolicy solver = - OdeSolverPolicy(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ r1, r2 }, parameters); + OdeSolverPolicy solver = OdeSolverPolicy( + micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ r1, r2 }, parameters); double temperature = 272.5; double pressure = 101253.3; @@ -747,7 +752,9 @@ void test_analytical_stiff_ternary_chemical_activation( .SetPhase(gas_phase); OdeSolverPolicy solver = OdeSolverPolicy( - micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ r1, r2, r3, r4, r5 }, parameters); + micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), + std::vector{ r1, r2, r3, r4, r5 }, + parameters); double temperature = 272.5; double pressure = 101253.3; @@ -857,8 +864,8 @@ void test_analytical_tunneling( .SetRateConstant(micm::TunnelingRateConstant({ .A_ = 1.2e-4, .B_ = 167, .C_ = 1.0e8 })) .SetPhase(gas_phase); - OdeSolverPolicy solver = - OdeSolverPolicy(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ r1, r2 }, parameters); + OdeSolverPolicy solver = OdeSolverPolicy( + micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ r1, r2 }, parameters); double temperature = 272.5; double pressure = 101253.3; @@ -982,7 +989,9 @@ void test_analytical_stiff_tunneling( .SetPhase(gas_phase); OdeSolverPolicy solver = OdeSolverPolicy( - micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ r1, r2, r3, r4, r5 }, parameters); + micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), + std::vector{ r1, r2, r3, r4, r5 }, + parameters); double temperature = 272.5; double pressure = 101253.3; @@ -1088,8 +1097,8 @@ void test_analytical_arrhenius( .SetRateConstant(micm::ArrheniusRateConstant({ .A_ = 1.2e-4, .B_ = 7, .C_ = 75, .D_ = 50, .E_ = 0.5 })) .SetPhase(gas_phase); - OdeSolverPolicy solver = - OdeSolverPolicy(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ r1, r2 }, parameters); + OdeSolverPolicy solver = OdeSolverPolicy( + micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ r1, r2 }, parameters); double temperature = 272.5; double pressure = 101253.3; @@ -1214,7 +1223,9 @@ void test_analytical_stiff_arrhenius( .SetPhase(gas_phase); OdeSolverPolicy solver = OdeSolverPolicy( - micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ r1, r2, r3, r4, r5 }, parameters); + micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), + std::vector{ r1, r2, r3, r4, r5 }, + parameters); double temperature = 272.5; double pressure = 101253.3; @@ -1330,8 +1341,8 @@ void test_analytical_branched( .n_ = 2 })) .SetPhase(gas_phase); - OdeSolverPolicy solver = - OdeSolverPolicy(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ r1, r2 }, parameters); + OdeSolverPolicy solver = OdeSolverPolicy( + micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ r1, r2 }, parameters); double temperature = 272.5; double pressure = 101253.3; @@ -1485,7 +1496,9 @@ void test_analytical_stiff_branched( .SetPhase(gas_phase); OdeSolverPolicy solver = OdeSolverPolicy( - micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ r1, r2, r3, r4, r5 }, parameters); + micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), + std::vector{ r1, r2, r3, r4, r5 }, + parameters); double temperature = 272.5; double pressure = 101253.3; diff --git a/test/integration/jit_terminator.cpp b/test/integration/jit_terminator.cpp index 254e6c5fa..180ca45ec 100644 --- a/test/integration/jit_terminator.cpp +++ b/test/integration/jit_terminator.cpp @@ -16,13 +16,11 @@ void RunTerminatorTest() solver_params.relative_tolerance_ = 1.0e-8; solver_params.max_number_of_steps_ = 100000; - TestTerminator< - micm::JitRosenbrockSolver< - MatrixPolicy, - SparseMatrixPolicy, - micm::JitLinearSolver, - micm::JitProcessSet>>( - number_of_grid_cells, solver_params); + TestTerminator, + micm::JitProcessSet>>(number_of_grid_cells, solver_params); } template diff --git a/test/integration/terminator.cpp b/test/integration/terminator.cpp index 0027e3463..3bd1913a4 100644 --- a/test/integration/terminator.cpp +++ b/test/integration/terminator.cpp @@ -19,10 +19,12 @@ void RunTerminatorTest(std::size_t number_of_grid_cells) solver_params.relative_tolerance_ = 1.0e-8; solver_params.max_number_of_steps_ = 100000; - TestTerminator>(number_of_grid_cells, solver_params); + TestTerminator>( + number_of_grid_cells, solver_params); solver_params.check_singularity_ = true; - TestTerminator< micm::RosenbrockSolver>(number_of_grid_cells, solver_params); + TestTerminator>( + number_of_grid_cells, solver_params); } TEST(RosenbrockSolver, Terminator) diff --git a/test/integration/terminator.hpp b/test/integration/terminator.hpp index 7a14babfb..4ff200b64 100644 --- a/test/integration/terminator.hpp +++ b/test/integration/terminator.hpp @@ -46,7 +46,9 @@ void TestTerminator( .SetRateConstant(micm::ArrheniusRateConstant({ .A_ = k2 })); auto solver = OdeSolverPolicy( - micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::vector{ toy_r1, toy_r2 }, parameters); + micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), + std::vector{ toy_r1, toy_r2 }, + parameters); auto state = solver.GetState(); auto get_double = std::bind(std::lognormal_distribution(-2.0, 2.0), std::default_random_engine()); diff --git a/test/regression/RosenbrockChapman/jit_util.hpp b/test/regression/RosenbrockChapman/jit_util.hpp index fe9e1caba..8b1ded4ac 100644 --- a/test/regression/RosenbrockChapman/jit_util.hpp +++ b/test/regression/RosenbrockChapman/jit_util.hpp @@ -2,12 +2,7 @@ #include -template< - template - class MatrixPolicy, - class SparseMatrixPolicy, - class LinearSolverPolicy, - class ProcessSetPolicy> +template class MatrixPolicy, class SparseMatrixPolicy, class LinearSolverPolicy, class ProcessSetPolicy> micm::JitRosenbrockSolver getTwoStageMultiCellJitChapmanSolver(const size_t number_of_grid_cells) { @@ -17,15 +12,11 @@ getTwoStageMultiCellJitChapmanSolver(const size_t number_of_grid_cells) auto options = micm::RosenbrockSolverParameters::TwoStageRosenbrockParameters(number_of_grid_cells); options.ignore_unused_species_ = true; - return micm::JitRosenbrockSolver(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::move(processes), options); + return micm::JitRosenbrockSolver( + micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::move(processes), options); } -template< - template - class MatrixPolicy, - class SparseMatrixPolicy, - class LinearSolverPolicy, - class ProcessSetPolicy> +template class MatrixPolicy, class SparseMatrixPolicy, class LinearSolverPolicy, class ProcessSetPolicy> micm::JitRosenbrockSolver getThreeStageMultiCellJitChapmanSolver(const size_t number_of_grid_cells) { @@ -35,15 +26,11 @@ getThreeStageMultiCellJitChapmanSolver(const size_t number_of_grid_cells) auto options = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters(number_of_grid_cells); options.ignore_unused_species_ = true; - return micm::JitRosenbrockSolver(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::move(processes), options); + return micm::JitRosenbrockSolver( + micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::move(processes), options); } -template< - template - class MatrixPolicy, - class SparseMatrixPolicy, - class LinearSolverPolicy, - class ProcessSetPolicy> +template class MatrixPolicy, class SparseMatrixPolicy, class LinearSolverPolicy, class ProcessSetPolicy> micm::JitRosenbrockSolver getFourStageMultiCellJitChapmanSolver(const size_t number_of_grid_cells) { @@ -53,15 +40,11 @@ getFourStageMultiCellJitChapmanSolver(const size_t number_of_grid_cells) auto options = micm::RosenbrockSolverParameters::FourStageRosenbrockParameters(number_of_grid_cells); options.ignore_unused_species_ = true; - return micm::JitRosenbrockSolver(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::move(processes), options); + return micm::JitRosenbrockSolver( + micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::move(processes), options); } -template< - template - class MatrixPolicy, - class SparseMatrixPolicy, - class LinearSolverPolicy, - class ProcessSetPolicy> +template class MatrixPolicy, class SparseMatrixPolicy, class LinearSolverPolicy, class ProcessSetPolicy> micm::JitRosenbrockSolver getFourStageDAMultiCellJitChapmanSolver(const size_t number_of_grid_cells) { @@ -71,15 +54,11 @@ getFourStageDAMultiCellJitChapmanSolver(const size_t number_of_grid_cells) auto options = micm::RosenbrockSolverParameters::FourStageDifferentialAlgebraicRosenbrockParameters(number_of_grid_cells); options.ignore_unused_species_ = true; - return micm::JitRosenbrockSolver(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::move(processes), options); + return micm::JitRosenbrockSolver( + micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::move(processes), options); } -template< - template - class MatrixPolicy, - class SparseMatrixPolicy, - class LinearSolverPolicy, - class ProcessSetPolicy> +template class MatrixPolicy, class SparseMatrixPolicy, class LinearSolverPolicy, class ProcessSetPolicy> micm::JitRosenbrockSolver getSixStageDAMultiCellJitChapmanSolver(const size_t number_of_grid_cells) { @@ -89,5 +68,6 @@ getSixStageDAMultiCellJitChapmanSolver(const size_t number_of_grid_cells) auto options = micm::RosenbrockSolverParameters::SixStageDifferentialAlgebraicRosenbrockParameters(number_of_grid_cells); options.ignore_unused_species_ = true; - return micm::JitRosenbrockSolver( micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::move(processes), options); + return micm::JitRosenbrockSolver( + micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }), std::move(processes), options); } \ No newline at end of file diff --git a/test/unit/process/test_cuda_process_set.cpp b/test/unit/process/test_cuda_process_set.cpp index ee0b84844..a9746b756 100644 --- a/test/unit/process/test_cuda_process_set.cpp +++ b/test/unit/process/test_cuda_process_set.cpp @@ -109,11 +109,7 @@ void testRandomSystemAddForcingTerms(std::size_t n_cells, std::size_t n_reaction } } -template< - class CPUMatrixPolicy, - class CPUSparseMatrixPolicy, - class GPUDenseMatrixPolicy, - class GPUSparseMatrixPolicy> +template void testRandomSystemSubtractJacobianTerms(std::size_t n_cells, std::size_t n_reactions, std::size_t n_species) { auto get_n_react = std::bind(std::uniform_int_distribution<>(0, 3), std::default_random_engine()); diff --git a/test/unit/process/test_jit_forcing_calculation.cpp b/test/unit/process/test_jit_forcing_calculation.cpp index 9ebe975ce..f5eb3ae7f 100644 --- a/test/unit/process/test_jit_forcing_calculation.cpp +++ b/test/unit/process/test_jit_forcing_calculation.cpp @@ -49,7 +49,7 @@ TEST(JitProcessSet, ForcingFunction) std::vector processes{ r1, r2, r3 }; - micm::JitProcessSet process_set{processes, state.variable_map_ }; + micm::JitProcessSet process_set{ processes, state.variable_map_ }; ForcingTestVectorMatrix rate_constants{ NUM_GRID_CELLS, 3 }; ForcingTestVectorMatrix forcing{ NUM_GRID_CELLS, 6, 0.0 }; diff --git a/test/unit/solver/test_backward_euler.cpp b/test/unit/solver/test_backward_euler.cpp index e85cd9113..5470d5660 100644 --- a/test/unit/solver/test_backward_euler.cpp +++ b/test/unit/solver/test_backward_euler.cpp @@ -5,7 +5,8 @@ #include #include -namespace { +namespace +{ auto a = micm::Species("A"); auto b = micm::Species("B"); auto c = micm::Species("C"); @@ -27,7 +28,7 @@ namespace { auto the_system = micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase }); std::vector reactions = { r1, r2 }; -} +} // namespace TEST(BackwardEuler, CanCallSolve) { @@ -35,11 +36,11 @@ TEST(BackwardEuler, CanCallSolve) params.absolute_tolerance_ = { 1e-6, 1e-6, 1e-6 }; auto be = micm::CpuSolverBuilder, micm::SparseMatrix>() - .SetSystem(the_system) - .SetReactions(reactions) - .SetNumberOfGridCells(1) - .SetSolverParameters(params) - .Build(); + .SetSystem(the_system) + .SetReactions(reactions) + .SetNumberOfGridCells(1) + .SetSolverParameters(params) + .Build(); double time_step = 1.0; auto state = be.GetState(); diff --git a/test/unit/solver/test_cuda_linear_solver.cpp b/test/unit/solver/test_cuda_linear_solver.cpp index eca79c22d..1fa3d2516 100644 --- a/test/unit/solver/test_cuda_linear_solver.cpp +++ b/test/unit/solver/test_cuda_linear_solver.cpp @@ -116,15 +116,18 @@ TEST(CudaLinearSolver, RandomMatrixVectorOrderingPolicy) testRandomMatrix>(1); testRandomMatrix>(20); testRandomMatrix>(300); - testRandomMatrix>(4000); + testRandomMatrix>( + 4000); } TEST(CudaLinearSolver, DiagonalMatrixVectorOrderingPolicy) { testDiagonalMatrix>(1); testDiagonalMatrix>(20); - testDiagonalMatrix>(300); - testDiagonalMatrix>(4000); + testDiagonalMatrix>( + 300); + testDiagonalMatrix>( + 4000); } TEST(CudaLinearSolver, RandomMatrixVectorOrderingForGPU) diff --git a/test/unit/solver/test_jit_linear_solver.cpp b/test/unit/solver/test_jit_linear_solver.cpp index 501152872..b0bc579d8 100644 --- a/test/unit/solver/test_jit_linear_solver.cpp +++ b/test/unit/solver/test_jit_linear_solver.cpp @@ -22,7 +22,10 @@ using Group4SparseVectorMatrix = micm::SparseMatrix>>(); + testDenseMatrix< + Group1VectorMatrix, + Group1SparseVectorMatrix, + micm::JitLinearSolver<1, Group1SparseVectorMatrix, micm::JitLuDecomposition<1>>>(); } TEST(JitLinearSolver, RandomMatrixVectorOrdering) diff --git a/test/unit/solver/test_linear_solver_policy.hpp b/test/unit/solver/test_linear_solver_policy.hpp index 9fe8f548d..488b2f27d 100644 --- a/test/unit/solver/test_linear_solver_policy.hpp +++ b/test/unit/solver/test_linear_solver_policy.hpp @@ -96,16 +96,16 @@ void testDenseMatrix() using FloatingPointType = typename MatrixPolicy::value_type; SparseMatrixPolicy A = SparseMatrixPolicy(SparseMatrixPolicy::Create(3) - .InitialValue(1.0e-30) - .WithElement(0, 0) - .WithElement(0, 1) - .WithElement(0, 2) - .WithElement(1, 0) - .WithElement(1, 1) - .WithElement(1, 2) - .WithElement(2, 0) - .WithElement(2, 1) - .WithElement(2, 2)); + .InitialValue(1.0e-30) + .WithElement(0, 0) + .WithElement(0, 1) + .WithElement(0, 2) + .WithElement(1, 0) + .WithElement(1, 1) + .WithElement(1, 2) + .WithElement(2, 0) + .WithElement(2, 1) + .WithElement(2, 2)); MatrixPolicy b(1, 3, 0.0); MatrixPolicy x(1, 3, 100.0); @@ -269,7 +269,7 @@ void testMarkowitzReordering() builder = builder.WithElement(i, j); SparseMatrixPolicy reordered_jac{ builder }; - auto orig_LU_calc = micm::LuDecomposition::Create< SparseMatrixPolicy>(orig_jac); + auto orig_LU_calc = micm::LuDecomposition::Create(orig_jac); auto reordered_LU_calc = micm::LuDecomposition::Create(reordered_jac); auto orig_LU = orig_LU_calc.template GetLUMatrices(orig_jac, 0.0); diff --git a/test/unit/solver/test_lu_decomposition_policy.hpp b/test/unit/solver/test_lu_decomposition_policy.hpp index eedc75c4b..f0c9fcb4b 100644 --- a/test/unit/solver/test_lu_decomposition_policy.hpp +++ b/test/unit/solver/test_lu_decomposition_policy.hpp @@ -77,16 +77,16 @@ template void testDenseMatrix() { SparseMatrixPolicy A = SparseMatrixPolicy(SparseMatrixPolicy::Create(3) - .InitialValue(1.0e-30) - .WithElement(0, 0) - .WithElement(0, 1) - .WithElement(0, 2) - .WithElement(1, 0) - .WithElement(1, 1) - .WithElement(1, 2) - .WithElement(2, 0) - .WithElement(2, 1) - .WithElement(2, 2)); + .InitialValue(1.0e-30) + .WithElement(0, 0) + .WithElement(0, 1) + .WithElement(0, 2) + .WithElement(1, 0) + .WithElement(1, 1) + .WithElement(1, 2) + .WithElement(2, 0) + .WithElement(2, 1) + .WithElement(2, 2)); A[0][0][0] = 2; A[0][0][1] = -1; @@ -108,12 +108,9 @@ void testDenseMatrix() template void testSingularMatrix() { - SparseMatrixPolicy A = SparseMatrixPolicy(SparseMatrixPolicy::Create(2) - .InitialValue(1.0e-30) - .WithElement(0, 0) - .WithElement(0, 1) - .WithElement(1, 0) - .WithElement(1, 1)); + SparseMatrixPolicy A = SparseMatrixPolicy( + SparseMatrixPolicy::Create(2).InitialValue(1.0e-30).WithElement(0, 0).WithElement(0, 1).WithElement(1, 0).WithElement( + 1, 1)); A[0][0][0] = 0; A[0][0][1] = 1;