diff --git a/include/micm/cuda/process/cuda_process_set.hpp b/include/micm/cuda/process/cuda_process_set.hpp index cee74f494..33aae6a18 100644 --- a/include/micm/cuda/process/cuda_process_set.hpp +++ b/include/micm/cuda/process/cuda_process_set.hpp @@ -52,25 +52,27 @@ namespace micm void SetJacobianFlatIds(const SparseMatrix& matrix); template - requires(CudaMatrix&& VectorizableDense) void AddForcingTerms( - const MatrixPolicy& rate_constants, - const MatrixPolicy& state_variables, - MatrixPolicy& forcing) const; + requires(CudaMatrix && VectorizableDense) + void AddForcingTerms(const MatrixPolicy& rate_constants, const MatrixPolicy& state_variables, MatrixPolicy& forcing) + const; template - requires(!CudaMatrix) void AddForcingTerms( - const MatrixPolicy& rate_constants, - const MatrixPolicy& state_variables, - MatrixPolicy& forcing) const; + requires(!CudaMatrix) + void AddForcingTerms(const MatrixPolicy& rate_constants, const MatrixPolicy& state_variables, MatrixPolicy& forcing) + const; template - requires( - CudaMatrix&& CudaMatrix&& VectorizableDense&& VectorizableSparse< - SparseMatrixPolicy>) void SubtractJacobianTerms(const MatrixPolicy& rate_constants, const MatrixPolicy& state_variables, SparseMatrixPolicy& jacobian) - const; + requires( + CudaMatrix && CudaMatrix && VectorizableDense && + VectorizableSparse) + void SubtractJacobianTerms( + const MatrixPolicy& rate_constants, + const MatrixPolicy& state_variables, + SparseMatrixPolicy& jacobian) const; template - requires(!CudaMatrix && !CudaMatrix) void SubtractJacobianTerms( + requires(!CudaMatrix && !CudaMatrix) + void SubtractJacobianTerms( const MatrixPolicy& rate_constants, const MatrixPolicy& state_variables, SparseMatrixPolicy& jacobian) const; @@ -119,7 +121,8 @@ namespace micm } template - requires(CudaMatrix&& VectorizableDense) inline void CudaProcessSet::AddForcingTerms( + requires(CudaMatrix && VectorizableDense) + inline void CudaProcessSet::AddForcingTerms( const MatrixPolicy& rate_constants, const MatrixPolicy& state_variables, MatrixPolicy& forcing) const @@ -130,12 +133,13 @@ namespace micm } template - requires(CudaMatrix&& CudaMatrix&& VectorizableDense&& - VectorizableSparse) inline void CudaProcessSet:: - SubtractJacobianTerms( - const MatrixPolicy& rate_constants, - const MatrixPolicy& state_variables, - SparseMatrixPolicy& jacobian) const + requires( + CudaMatrix && CudaMatrix && VectorizableDense && + VectorizableSparse) + inline void CudaProcessSet::SubtractJacobianTerms( + const MatrixPolicy& rate_constants, + const MatrixPolicy& state_variables, + SparseMatrixPolicy& jacobian) const { auto jacobian_param = jacobian.AsDeviceParam(); // we need to update jacobian so it can't be constant and must be an lvalue diff --git a/include/micm/cuda/solver/cuda_linear_solver.hpp b/include/micm/cuda/solver/cuda_linear_solver.hpp index a84083f1e..657f6ff2e 100644 --- a/include/micm/cuda/solver/cuda_linear_solver.hpp +++ b/include/micm/cuda/solver/cuda_linear_solver.hpp @@ -79,9 +79,10 @@ namespace micm }; template - requires( - CudaMatrix&& CudaMatrix&& VectorizableDense&& VectorizableSparse< - SparseMatrixPolicy>) void Solve(MatrixPolicy& x, const SparseMatrixPolicy& L, const SparseMatrixPolicy& U) const + requires( + CudaMatrix && CudaMatrix && VectorizableDense && + VectorizableSparse) + void Solve(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 micm::cuda::SolveKernelDriver(x_param, L.AsDeviceParam(), U.AsDeviceParam(), this->devstruct_); diff --git a/include/micm/cuda/solver/cuda_lu_decomposition_doolittle.hpp b/include/micm/cuda/solver/cuda_lu_decomposition_doolittle.hpp index f3f7edff9..978bbe647 100644 --- a/include/micm/cuda/solver/cuda_lu_decomposition_doolittle.hpp +++ b/include/micm/cuda/solver/cuda_lu_decomposition_doolittle.hpp @@ -85,10 +85,15 @@ namespace micm /// @brief Create an LU decomposition algorithm for a given sparse matrix policy /// @param matrix Sparse matrix template - requires(SparseMatrixConcept) static CudaLuDecompositionDoolittle Create(const SparseMatrixPolicy& matrix) + requires(SparseMatrixConcept) + static CudaLuDecompositionDoolittle Create(const SparseMatrixPolicy& matrix) { - static_assert(std::is_same_v, "SparseMatrixPolicy must be the same as LMatrixPolicy for CUDA LU decomposition"); - static_assert(std::is_same_v, "SparseMatrixPolicy must be the same as UMatrixPolicy for CUDA LU decomposition"); + static_assert( + std::is_same_v, + "SparseMatrixPolicy must be the same as LMatrixPolicy for CUDA LU decomposition"); + static_assert( + std::is_same_v, + "SparseMatrixPolicy must be the same as UMatrixPolicy for CUDA LU decomposition"); CudaLuDecompositionDoolittle lu_decomp(matrix); return lu_decomp; } @@ -98,15 +103,13 @@ namespace micm /// @param L is the lower triangular matrix created by decomposition /// @param U is the upper triangular matrix created by decomposition template - requires(CudaMatrix&& VectorizableSparse) void Decompose( - const SparseMatrixPolicy& A, - auto& L, - auto& U) const; + requires(CudaMatrix && VectorizableSparse) + void Decompose(const SparseMatrixPolicy& A, auto& L, auto& U) const; }; template - requires(CudaMatrix&& VectorizableSparse) void CudaLuDecompositionDoolittle:: - Decompose(const SparseMatrixPolicy& A, auto& L, auto& U) const + requires(CudaMatrix && VectorizableSparse) + void CudaLuDecompositionDoolittle::Decompose(const SparseMatrixPolicy& A, auto& L, auto& 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/cuda/solver/cuda_rosenbrock.hpp b/include/micm/cuda/solver/cuda_rosenbrock.hpp index 2044f97c9..836533d7c 100644 --- a/include/micm/cuda/solver/cuda_rosenbrock.hpp +++ b/include/micm/cuda/solver/cuda_rosenbrock.hpp @@ -101,7 +101,7 @@ namespace micm /// @param alpha template void AlphaMinusJacobian(SparseMatrixPolicy& jacobian, const double& alpha) const - requires(CudaMatrix&& VectorizableSparse) + requires(CudaMatrix && VectorizableSparse) { auto jacobian_param = jacobian.AsDeviceParam(); // we need to update jacobian so it can't be constant and must be an lvalue @@ -116,7 +116,8 @@ namespace micm /// @return The scaled norm of the errors template double NormalizedError(const DenseMatrixPolicy& y_old, const DenseMatrixPolicy& y_new, const DenseMatrixPolicy& errors) - const requires(CudaMatrix&& VectorizableDense) + const + requires(CudaMatrix && VectorizableDense) { return micm::cuda::NormalizedErrorDriver( y_old.AsDeviceParam(), y_new.AsDeviceParam(), errors.AsDeviceParam(), this->parameters_, this->devstruct_); diff --git a/include/micm/cuda/solver/cuda_state.hpp b/include/micm/cuda/solver/cuda_state.hpp index a0b0fc069..a20ad1a29 100644 --- a/include/micm/cuda/solver/cuda_state.hpp +++ b/include/micm/cuda/solver/cuda_state.hpp @@ -25,14 +25,16 @@ namespace micm : State(parameters){}; /// @brief Copy input variables to the device - void SyncInputsToDevice() requires(CudaMatrix&& VectorizableDense) + void SyncInputsToDevice() + requires(CudaMatrix && VectorizableDense) { this->variables_.CopyToDevice(); this->rate_constants_.CopyToDevice(); } /// @brief Copy output variables to the host - void SyncOutputsToHost() requires(CudaMatrix&& VectorizableDense) + void SyncOutputsToHost() + requires(CudaMatrix && VectorizableDense) { this->variables_.CopyToHost(); } diff --git a/include/micm/cuda/util/cuda_dense_matrix.hpp b/include/micm/cuda/util/cuda_dense_matrix.hpp index 13b3552da..d0994c3c9 100644 --- a/include/micm/cuda/util/cuda_dense_matrix.hpp +++ b/include/micm/cuda/util/cuda_dense_matrix.hpp @@ -42,17 +42,10 @@ namespace micm /// Concept for Cuda Matrix template - concept CudaMatrix = requires(MatrixType t) - { - { - t.CopyToDevice() - } -> std::same_as; - { - t.CopyToHost() - } -> std::same_as; - { - t.AsDeviceParam() - } -> std::same_as; + concept CudaMatrix = requires(MatrixType t) { + { t.CopyToDevice() } -> std::same_as; + { t.CopyToHost() } -> std::same_as; + { t.AsDeviceParam() } -> std::same_as; }; template @@ -192,6 +185,22 @@ namespace micm "CUBLAS Daxpy operation failed..."); } + /// @brief For each element of the VectorMatrix, perform y = max(y, x), where x is a scalar constant + /// @param x The scalar constant to compare against + void Max(const T x) + { + static_assert(std::is_same_v); + CHECK_CUDA_ERROR(micm::cuda::MatrixMax(this->param_, x), "CudaMatrixMax"); + } + + /// @brief For each element of the VectorMatrix, perform y = min(y, x), where x is a scalar constant + /// @param x The scalar constant to compare against + void Min(const T x) + { + static_assert(std::is_same_v); + CHECK_CUDA_ERROR(micm::cuda::MatrixMin(this->param_, x), "CudaMatrixMin"); + } + // Copy the device data from the other Cuda dense matrix into this one void Copy(const CudaDenseMatrix& other) { diff --git a/include/micm/cuda/util/cuda_matrix.cuh b/include/micm/cuda/util/cuda_matrix.cuh index 61bfff6e5..713f5844c 100644 --- a/include/micm/cuda/util/cuda_matrix.cuh +++ b/include/micm/cuda/util/cuda_matrix.cuh @@ -24,6 +24,18 @@ namespace micm /// @returns Error code from free-ing data on device, if any cudaError_t FreeVector(CudaMatrixParam& vectorMatrix); + /// @brief Sets all elements in a CUDA matrix to their current value or a specified value, whichever is greater + /// @param vectorMatrix Struct containing allocated device memory + /// @param val Value to compare with each element in the matrix + template + cudaError_t MatrixMax(CudaMatrixParam& vectorMatrix, T val); + + /// @brief Sets all elements in a CUDA matrix to their current value or a specified value, whichever is lesser + /// @param vectorMatrix Struct containing allocated device memory + /// @param val Value to compare with each element in the matrix + template + cudaError_t MatrixMin(CudaMatrixParam& vectorMatrix, T val); + /// @brief Copies data from the host to the device /// @param vectorMatrix Struct containing allocated device memory /// @param h_data Host data to copy from diff --git a/include/micm/jit/solver/jit_linear_solver.inl b/include/micm/jit/solver/jit_linear_solver.inl index b678aeac1..f4e401f0f 100644 --- a/include/micm/jit/solver/jit_linear_solver.inl +++ b/include/micm/jit/solver/jit_linear_solver.inl @@ -14,8 +14,8 @@ namespace micm } template - inline JitLinearSolver - &JitLinearSolver::operator=(JitLinearSolver &&other) + inline JitLinearSolver & + JitLinearSolver::operator=(JitLinearSolver &&other) { LinearSolver::operator=(std::move(other)); solve_function_resource_tracker_ = std::move(other.solve_function_resource_tracker_); diff --git a/include/micm/jit/solver/jit_lu_decomposition_doolittle.hpp b/include/micm/jit/solver/jit_lu_decomposition_doolittle.hpp index 3045b7789..aec7421dd 100644 --- a/include/micm/jit/solver/jit_lu_decomposition_doolittle.hpp +++ b/include/micm/jit/solver/jit_lu_decomposition_doolittle.hpp @@ -38,14 +38,15 @@ namespace micm /// @brief Create an LU decomposition algorithm for a given sparse matrix policy /// @param matrix Sparse matrix template - requires(SparseMatrixConcept) static JitLuDecompositionDoolittle Create(const SparseMatrixPolicy& matrix); + requires(SparseMatrixConcept) + static JitLuDecompositionDoolittle Create(const SparseMatrixPolicy &matrix); /// @brief Create sparse L and U matrices for a given A matrix /// @param A Sparse matrix that will be decomposed /// @param lower The lower triangular matrix created by decomposition /// @param upper The upper triangular matrix created by decomposition template - void Decompose(const SparseMatrixPolicy &A, auto& lower, auto& upper) const; + void Decompose(const SparseMatrixPolicy &A, auto &lower, auto &upper) const; private: /// @brief Generates a function to perform the LU decomposition for a specific matrix sparsity structure diff --git a/include/micm/jit/solver/jit_lu_decomposition_doolittle.inl b/include/micm/jit/solver/jit_lu_decomposition_doolittle.inl index 95cc27590..4a208d77e 100644 --- a/include/micm/jit/solver/jit_lu_decomposition_doolittle.inl +++ b/include/micm/jit/solver/jit_lu_decomposition_doolittle.inl @@ -38,7 +38,8 @@ namespace micm std::to_string(L); throw std::system_error(make_error_code(MicmJitErrc::InvalidMatrix), msg); } - Initialize(matrix, typename SparseMatrixPolicy::value_type()); + Initialize( + matrix, typename SparseMatrixPolicy::value_type()); GenerateDecomposeFunction(); } @@ -54,12 +55,15 @@ namespace micm template template - requires(SparseMatrixConcept) - inline JitLuDecompositionDoolittle JitLuDecompositionDoolittle::Create( - const SparseMatrixPolicy& matrix) + requires(SparseMatrixConcept) + inline JitLuDecompositionDoolittle JitLuDecompositionDoolittle::Create(const SparseMatrixPolicy &matrix) { - static_assert(std::is_same_v, "SparseMatrixPolicy must be the same as LMatrixPolicy for JIT LU decomposition"); - static_assert(std::is_same_v, "SparseMatrixPolicy must be the same as UMatrixPolicy for JIT LU decomposition"); + static_assert( + std::is_same_v, + "SparseMatrixPolicy must be the same as LMatrixPolicy for JIT LU decomposition"); + static_assert( + std::is_same_v, + "SparseMatrixPolicy must be the same as UMatrixPolicy for JIT LU decomposition"); JitLuDecompositionDoolittle lu_decomp(matrix); return lu_decomp; } @@ -220,10 +224,7 @@ namespace micm template template - void JitLuDecompositionDoolittle::Decompose( - const SparseMatrixPolicy &A, - auto& lower, - auto& upper) const + void JitLuDecompositionDoolittle::Decompose(const SparseMatrixPolicy &A, auto &lower, auto &upper) const { decompose_function_(A.AsVector().data(), lower.AsVector().data(), upper.AsVector().data()); for (size_t block = 0; block < A.NumberOfBlocks(); ++block) diff --git a/include/micm/process/process.hpp b/include/micm/process/process.hpp index 9d0868ef1..dacdec3be 100644 --- a/include/micm/process/process.hpp +++ b/include/micm/process/process.hpp @@ -84,12 +84,24 @@ namespace micm /// @brief Recalculate the rate constants for each process for the current state /// @param processes The set of processes being solved /// @param state The solver state to update - template - requires(!VectorizableDense) static void CalculateRateConstants( + template< + class DenseMatrixPolicy, + class SparseMatrixPolicy, + class LuDecompositionPolicy, + class LMatrixPolicy, + class UMatrixPolicy> + requires(!VectorizableDense) + static void CalculateRateConstants( const std::vector& processes, State& state); - template - requires(VectorizableDense) static void CalculateRateConstants( + template< + class DenseMatrixPolicy, + class SparseMatrixPolicy, + class LuDecompositionPolicy, + class LMatrixPolicy, + class UMatrixPolicy> + requires(VectorizableDense) + static void CalculateRateConstants( const std::vector& processes, State& state); @@ -147,8 +159,14 @@ namespace micm ProcessBuilder& SetPhase(const Phase& phase); }; - template - requires(!VectorizableDense) void Process::CalculateRateConstants( + template< + class DenseMatrixPolicy, + class SparseMatrixPolicy, + class LuDecompositionPolicy, + class LMatrixPolicy, + class UMatrixPolicy> + requires(!VectorizableDense) + void Process::CalculateRateConstants( const std::vector& processes, State& state) { @@ -172,8 +190,14 @@ namespace micm } } - template - requires(VectorizableDense) void Process::CalculateRateConstants( + template< + class DenseMatrixPolicy, + class SparseMatrixPolicy, + class LuDecompositionPolicy, + class LMatrixPolicy, + class UMatrixPolicy> + requires(VectorizableDense) + void Process::CalculateRateConstants( const std::vector& processes, State& state) { diff --git a/include/micm/process/process_set.hpp b/include/micm/process/process_set.hpp index 9cb20c90f..44985d13e 100644 --- a/include/micm/process/process_set.hpp +++ b/include/micm/process/process_set.hpp @@ -92,12 +92,13 @@ namespace micm /// @param state_variables Current state variable values (grid cell, state variable) /// @param forcing Forcing terms for each state variable (grid cell, state variable) template - requires(!VectorizableDense) void AddForcingTerms( + requires(!VectorizableDense) + void AddForcingTerms( const DenseMatrixPolicy& rate_constants, const DenseMatrixPolicy& state_variables, DenseMatrixPolicy& forcing) const; template - requires VectorizableDense + requires VectorizableDense void AddForcingTerms( const DenseMatrixPolicy& rate_constants, const DenseMatrixPolicy& state_variables, @@ -108,13 +109,15 @@ namespace micm /// @param state_variables Current state variable values (grid cell, state variable) /// @param jacobian Jacobian matrix for the system (grid cell, dependent variable, independent variable) template - requires(!VectorizableDense || !VectorizableSparse) void SubtractJacobianTerms( + requires(!VectorizableDense || !VectorizableSparse) + void SubtractJacobianTerms( const DenseMatrixPolicy& rate_constants, const DenseMatrixPolicy& state_variables, SparseMatrixPolicy& jacobian) const; // template class MatrixPolicy, template class SparseMatrixPolicy> template - requires(VectorizableDense&& VectorizableSparse) void SubtractJacobianTerms( + requires(VectorizableDense && VectorizableSparse) + void SubtractJacobianTerms( const DenseMatrixPolicy& rate_constants, const DenseMatrixPolicy& state_variables, SparseMatrixPolicy& jacobian) const; @@ -217,7 +220,8 @@ namespace micm } template - requires(!VectorizableDense) inline void ProcessSet::AddForcingTerms( + requires(!VectorizableDense) + inline void ProcessSet::AddForcingTerms( const DenseMatrixPolicy& rate_constants, const DenseMatrixPolicy& state_variables, DenseMatrixPolicy& forcing) const @@ -261,7 +265,7 @@ namespace micm }; template - requires VectorizableDense + requires VectorizableDense inline void ProcessSet::AddForcingTerms( const DenseMatrixPolicy& rate_constants, const DenseMatrixPolicy& state_variables, @@ -306,11 +310,11 @@ namespace micm // Forming the Jacobian matrix "J" and returning "-J" to be consistent with the CUDA implementation template - requires(!VectorizableDense || !VectorizableSparse) inline void ProcessSet:: - SubtractJacobianTerms( - const DenseMatrixPolicy& rate_constants, - const DenseMatrixPolicy& state_variables, - SparseMatrixPolicy& jacobian) const + requires(!VectorizableDense || !VectorizableSparse) + inline void ProcessSet::SubtractJacobianTerms( + const DenseMatrixPolicy& rate_constants, + const DenseMatrixPolicy& state_variables, + SparseMatrixPolicy& jacobian) const { MICM_PROFILE_FUNCTION(); @@ -361,11 +365,11 @@ namespace micm // Forming the Jacobian matrix "J" and returning "-J" to be consistent with the CUDA implementation template - requires(VectorizableDense&& VectorizableSparse) inline void ProcessSet:: - SubtractJacobianTerms( - const DenseMatrixPolicy& rate_constants, - const DenseMatrixPolicy& state_variables, - SparseMatrixPolicy& jacobian) const + requires(VectorizableDense && VectorizableSparse) + inline void ProcessSet::SubtractJacobianTerms( + const DenseMatrixPolicy& rate_constants, + const DenseMatrixPolicy& state_variables, + SparseMatrixPolicy& jacobian) const { MICM_PROFILE_FUNCTION(); diff --git a/include/micm/solver/backward_euler.hpp b/include/micm/solver/backward_euler.hpp index 808ce8d05..30424a6a3 100644 --- a/include/micm/solver/backward_euler.hpp +++ b/include/micm/solver/backward_euler.hpp @@ -80,12 +80,14 @@ namespace micm static bool IsConverged( const BackwardEulerSolverParameters& parameters, const DenseMatrixPolicy& residual, - const DenseMatrixPolicy& state) requires(!VectorizableDense); + const DenseMatrixPolicy& state) + requires(!VectorizableDense); template static bool IsConverged( const BackwardEulerSolverParameters& parameters, const DenseMatrixPolicy& residual, - const DenseMatrixPolicy& state) requires(VectorizableDense); + const DenseMatrixPolicy& state) + requires(VectorizableDense); }; } // namespace micm diff --git a/include/micm/solver/backward_euler.inl b/include/micm/solver/backward_euler.inl index 620b82efe..b810a7c8d 100644 --- a/include/micm/solver/backward_euler.inl +++ b/include/micm/solver/backward_euler.inl @@ -195,7 +195,8 @@ namespace micm inline bool BackwardEuler::IsConverged( const BackwardEulerSolverParameters& parameters, const DenseMatrixPolicy& residual, - const DenseMatrixPolicy& state) requires(!VectorizableDense) + const DenseMatrixPolicy& state) + requires(!VectorizableDense) { double small = parameters.small_; double rel_tol = parameters.relative_tolerance_; @@ -221,7 +222,8 @@ namespace micm inline bool BackwardEuler::IsConverged( const BackwardEulerSolverParameters& parameters, const DenseMatrixPolicy& residual, - const DenseMatrixPolicy& state) requires(VectorizableDense) + const DenseMatrixPolicy& state) + requires(VectorizableDense) { double small = parameters.small_; double rel_tol = parameters.relative_tolerance_; diff --git a/include/micm/solver/linear_solver.hpp b/include/micm/solver/linear_solver.hpp index 76396cfe6..34631fd9a 100644 --- a/include/micm/solver/linear_solver.hpp +++ b/include/micm/solver/linear_solver.hpp @@ -41,7 +41,11 @@ namespace micm /// @brief A general-use block-diagonal sparse-matrix linear solver /// /// The sparsity pattern of each block in the block diagonal matrix is the same. - template + template< + class SparseMatrixPolicy, + class LuDecompositionPolicy = LuDecomposition, + class LMatrixPolicy = SparseMatrixPolicy, + class UMatrixPolicy = SparseMatrixPolicy> class LinearSolver { protected: @@ -101,15 +105,11 @@ namespace micm /// @brief Solve for x in Ax = b. x should be a copy of b and after Solve finishes x will contain the result template - requires(!VectorizableDense || !VectorizableSparse) void Solve( - MatrixPolicy& x, - const LMatrixPolicy& lower_matrix, - const UMatrixPolicy& upper_matrix) const; + requires(!VectorizableDense || !VectorizableSparse) + void Solve(MatrixPolicy& x, const LMatrixPolicy& lower_matrix, const UMatrixPolicy& upper_matrix) const; template - requires(VectorizableDense&& VectorizableSparse) void Solve( - MatrixPolicy& x, - const LMatrixPolicy& lower_matrix, - const UMatrixPolicy& upper_matrix) const; + requires(VectorizableDense && VectorizableSparse) + void Solve(MatrixPolicy& x, const LMatrixPolicy& lower_matrix, const UMatrixPolicy& upper_matrix) const; }; } // namespace micm diff --git a/include/micm/solver/linear_solver.inl b/include/micm/solver/linear_solver.inl index 2efde35fd..ee774c1b7 100644 --- a/include/micm/solver/linear_solver.inl +++ b/include/micm/solver/linear_solver.inl @@ -118,10 +118,11 @@ namespace micm template template - requires( - !VectorizableDense || - !VectorizableSparse) inline void LinearSolver:: - Solve(MatrixPolicy& x, const LMatrixPolicy& lower_matrix, const UMatrixPolicy& upper_matrix) const + requires(!VectorizableDense || !VectorizableSparse) + inline void LinearSolver::Solve( + MatrixPolicy& x, + const LMatrixPolicy& lower_matrix, + const UMatrixPolicy& upper_matrix) const { MICM_PROFILE_FUNCTION(); @@ -173,9 +174,11 @@ namespace micm template template - requires(VectorizableDense&& - VectorizableSparse) inline void LinearSolver:: - Solve(MatrixPolicy& x, const LMatrixPolicy& lower_matrix, const UMatrixPolicy& upper_matrix) const + requires(VectorizableDense && VectorizableSparse) + inline void LinearSolver::Solve( + MatrixPolicy& x, + const LMatrixPolicy& lower_matrix, + const UMatrixPolicy& upper_matrix) const { MICM_PROFILE_FUNCTION(); constexpr std::size_t n_cells = MatrixPolicy::GroupVectorSize(); @@ -183,10 +186,8 @@ namespace micm for (std::size_t i_group = 0; i_group < x.NumberOfGroups(); ++i_group) { 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()); - auto U_group = - std::next(upper_matrix.AsVector().begin(), i_group * upper_matrix.GroupSize()); + auto L_group = std::next(lower_matrix.AsVector().begin(), i_group * lower_matrix.GroupSize()); + auto U_group = std::next(upper_matrix.AsVector().begin(), i_group * upper_matrix.GroupSize()); // Forward Substitution { auto y_elem = x_group; diff --git a/include/micm/solver/lu_decomposition_doolittle.hpp b/include/micm/solver/lu_decomposition_doolittle.hpp index ef4d98b45..4e6e482d8 100644 --- a/include/micm/solver/lu_decomposition_doolittle.hpp +++ b/include/micm/solver/lu_decomposition_doolittle.hpp @@ -90,20 +90,23 @@ namespace micm /// @brief Construct an LU decomposition algorithm for a given sparse matrix /// @param matrix Sparse matrix template - requires(SparseMatrixConcept) LuDecompositionDoolittle(const SparseMatrixPolicy& matrix); + requires(SparseMatrixConcept) + LuDecompositionDoolittle(const SparseMatrixPolicy& matrix); ~LuDecompositionDoolittle() = default; /// @brief Create an LU decomposition algorithm for a given sparse matrix policy /// @param matrix Sparse matrix template - requires(SparseMatrixConcept) static LuDecompositionDoolittle Create(const SparseMatrixPolicy& matrix); + requires(SparseMatrixConcept) + static LuDecompositionDoolittle Create(const SparseMatrixPolicy& matrix); /// @brief Create sparse L and U matrices for a given A matrix /// @param A Sparse matrix that will be decomposed /// @return L and U Sparse matrices template - requires(SparseMatrixConcept) static std::pair GetLUMatrices( + requires(SparseMatrixConcept) + static std::pair GetLUMatrices( const SparseMatrixPolicy& A, typename SparseMatrixPolicy::value_type initial_value); @@ -112,21 +115,18 @@ namespace micm /// @param L The lower triangular matrix created by decomposition /// @param U The upper triangular matrix created by decomposition template - requires(!VectorizableSparse) void Decompose( - const SparseMatrixPolicy& A, - auto& L, - auto& U) const; + requires(!VectorizableSparse) + void Decompose(const SparseMatrixPolicy& A, auto& L, auto& U) const; template - requires(VectorizableSparse) void Decompose( - const SparseMatrixPolicy& A, - auto& L, - auto& U) const; + requires(VectorizableSparse) + void Decompose(const SparseMatrixPolicy& A, auto& L, auto& U) const; protected: /// @brief Initialize arrays for the LU decomposition /// @param A Sparse matrix to decompose template - requires(SparseMatrixConcept) void Initialize(const SparseMatrixPolicy& matrix, auto initial_value); + requires(SparseMatrixConcept) + void Initialize(const SparseMatrixPolicy& matrix, auto initial_value); }; } // namespace micm diff --git a/include/micm/solver/lu_decomposition_doolittle.inl b/include/micm/solver/lu_decomposition_doolittle.inl index efd5147f4..512a1c68a 100644 --- a/include/micm/solver/lu_decomposition_doolittle.inl +++ b/include/micm/solver/lu_decomposition_doolittle.inl @@ -9,24 +9,25 @@ namespace micm } template - requires(SparseMatrixConcept) inline LuDecompositionDoolittle::LuDecompositionDoolittle(const SparseMatrixPolicy& matrix) + requires(SparseMatrixConcept) + inline LuDecompositionDoolittle::LuDecompositionDoolittle(const SparseMatrixPolicy& matrix) { Initialize(matrix, typename SparseMatrixPolicy::value_type()); } template - requires(SparseMatrixConcept) inline LuDecompositionDoolittle LuDecompositionDoolittle::Create( - const SparseMatrixPolicy& matrix) + requires(SparseMatrixConcept) + inline LuDecompositionDoolittle LuDecompositionDoolittle::Create(const SparseMatrixPolicy& matrix) { LuDecompositionDoolittle lu_decomp{}; - lu_decomp.Initialize(matrix, typename SparseMatrixPolicy::value_type()); + lu_decomp.Initialize( + matrix, typename SparseMatrixPolicy::value_type()); return lu_decomp; } template - requires(SparseMatrixConcept) inline void LuDecompositionDoolittle::Initialize( - const SparseMatrixPolicy& matrix, - auto initial_value) + requires(SparseMatrixConcept) + inline void LuDecompositionDoolittle::Initialize(const SparseMatrixPolicy& matrix, auto initial_value) { MICM_PROFILE_FUNCTION(); @@ -93,9 +94,10 @@ namespace micm } template - requires( - SparseMatrixConcept) inline std::pair LuDecompositionDoolittle:: - GetLUMatrices(const SparseMatrixPolicy& A, typename SparseMatrixPolicy::value_type initial_value) + requires(SparseMatrixConcept) + inline std::pair LuDecompositionDoolittle::GetLUMatrices( + const SparseMatrixPolicy& A, + typename SparseMatrixPolicy::value_type initial_value) { MICM_PROFILE_FUNCTION(); @@ -153,10 +155,8 @@ namespace micm } template - requires(!VectorizableSparse) inline void LuDecompositionDoolittle::Decompose( - const SparseMatrixPolicy& A, - auto& L, - auto& U) const + requires(!VectorizableSparse) + inline void LuDecompositionDoolittle::Decompose(const SparseMatrixPolicy& A, auto& L, auto& U) const { MICM_PROFILE_FUNCTION(); @@ -221,10 +221,8 @@ namespace micm } template - requires(VectorizableSparse) inline void LuDecompositionDoolittle::Decompose( - const SparseMatrixPolicy& A, - auto& L, - auto& U) const + requires(VectorizableSparse) + inline void LuDecompositionDoolittle::Decompose(const SparseMatrixPolicy& A, auto& L, auto& U) const { MICM_PROFILE_FUNCTION(); diff --git a/include/micm/solver/lu_decomposition_mozart.hpp b/include/micm/solver/lu_decomposition_mozart.hpp index 550f4e9da..39aad214c 100644 --- a/include/micm/solver/lu_decomposition_mozart.hpp +++ b/include/micm/solver/lu_decomposition_mozart.hpp @@ -96,20 +96,23 @@ namespace micm /// @brief Construct an LU decomposition algorithm for a given sparse matrix /// @param matrix Sparse matrix template - requires(SparseMatrixConcept) LuDecompositionMozart(const SparseMatrixPolicy& matrix); + requires(SparseMatrixConcept) + LuDecompositionMozart(const SparseMatrixPolicy& matrix); ~LuDecompositionMozart() = default; /// @brief Create an LU decomposition algorithm for a given sparse matrix policy /// @param matrix Sparse matrix template - requires(SparseMatrixConcept) static LuDecompositionMozart Create(const SparseMatrixPolicy& matrix); + requires(SparseMatrixConcept) + static LuDecompositionMozart Create(const SparseMatrixPolicy& matrix); /// @brief Create sparse L and U matrices for a given A matrix /// @param A Sparse matrix that will be decomposed /// @return L and U Sparse matrices template - requires(SparseMatrixConcept) static std::pair GetLUMatrices( + requires(SparseMatrixConcept) + static std::pair GetLUMatrices( const SparseMatrixPolicy& A, typename SparseMatrixPolicy::value_type initial_value); @@ -118,21 +121,18 @@ namespace micm /// @param L The lower triangular matrix created by decomposition /// @param U The upper triangular matrix created by decomposition template - requires(!VectorizableSparse) void Decompose( - const SparseMatrixPolicy& A, - auto& L, - auto& U) const; + requires(!VectorizableSparse) + void Decompose(const SparseMatrixPolicy& A, auto& L, auto& U) const; template - requires(VectorizableSparse) void Decompose( - const SparseMatrixPolicy& A, - auto& L, - auto& U) const; + requires(VectorizableSparse) + void Decompose(const SparseMatrixPolicy& A, auto& L, auto& U) const; private: /// @brief Initialize arrays for the LU decomposition /// @param A Sparse matrix to decompose template - requires(SparseMatrixConcept) void Initialize(const SparseMatrixPolicy& matrix, auto initial_value); + requires(SparseMatrixConcept) + void Initialize(const SparseMatrixPolicy& matrix, auto initial_value); }; } // namespace micm diff --git a/include/micm/solver/lu_decomposition_mozart.inl b/include/micm/solver/lu_decomposition_mozart.inl index c8bcd7f1e..c222bf682 100644 --- a/include/micm/solver/lu_decomposition_mozart.inl +++ b/include/micm/solver/lu_decomposition_mozart.inl @@ -9,24 +9,25 @@ namespace micm } template - requires(SparseMatrixConcept) inline LuDecompositionMozart::LuDecompositionMozart(const SparseMatrixPolicy& matrix) + requires(SparseMatrixConcept) + inline LuDecompositionMozart::LuDecompositionMozart(const SparseMatrixPolicy& matrix) { Initialize(matrix, typename SparseMatrixPolicy::value_type()); } template - requires(SparseMatrixConcept) inline LuDecompositionMozart LuDecompositionMozart::Create( - const SparseMatrixPolicy& matrix) + requires(SparseMatrixConcept) + inline LuDecompositionMozart LuDecompositionMozart::Create(const SparseMatrixPolicy& matrix) { LuDecompositionMozart lu_decomp{}; - lu_decomp.Initialize(matrix, typename SparseMatrixPolicy::value_type()); + lu_decomp.Initialize( + matrix, typename SparseMatrixPolicy::value_type()); return lu_decomp; } template - requires(SparseMatrixConcept) inline void LuDecompositionMozart::Initialize( - const SparseMatrixPolicy& matrix, - auto initial_value) + requires(SparseMatrixConcept) + inline void LuDecompositionMozart::Initialize(const SparseMatrixPolicy& matrix, auto initial_value) { MICM_PROFILE_FUNCTION(); @@ -111,9 +112,10 @@ namespace micm } template - requires( - SparseMatrixConcept) inline std::pair LuDecompositionMozart:: - GetLUMatrices(const SparseMatrixPolicy& A, typename SparseMatrixPolicy::value_type initial_value) + requires(SparseMatrixConcept) + inline std::pair LuDecompositionMozart::GetLUMatrices( + const SparseMatrixPolicy& A, + typename SparseMatrixPolicy::value_type initial_value) { MICM_PROFILE_FUNCTION(); @@ -163,10 +165,8 @@ namespace micm } template - requires(!VectorizableSparse) inline void LuDecompositionMozart::Decompose( - const SparseMatrixPolicy& A, - auto& L, - auto& U) const + requires(!VectorizableSparse) + inline void LuDecompositionMozart::Decompose(const SparseMatrixPolicy& A, auto& L, auto& U) const { MICM_PROFILE_FUNCTION(); const std::size_t n = A.NumRows(); @@ -184,7 +184,7 @@ namespace micm auto nujk_nljk_uik = nujk_nljk_uik_.begin(); auto ujk_lji = ujk_lji_.begin(); auto ljk_lji = ljk_lji_.begin(); - + for (auto& lii_nuji_nlji : lii_nuji_nlji_) { for (std::size_t i = 0; i < std::get<1>(lii_nuji_nlji); ++i) @@ -232,10 +232,8 @@ namespace micm } template - requires(VectorizableSparse) inline void LuDecompositionMozart::Decompose( - const SparseMatrixPolicy& A, - auto& L, - auto& U) const + requires(VectorizableSparse) + inline void LuDecompositionMozart::Decompose(const SparseMatrixPolicy& A, auto& L, auto& U) const { MICM_PROFILE_FUNCTION(); @@ -266,24 +264,24 @@ namespace micm for (std::size_t i = 0; i < std::get<1>(lii_nuji_nlji); ++i) { for (std::size_t i_cell = 0; i_cell < n_cells; ++i_cell) - U_vector[uji_aji->first+i_cell] = A_vector[uji_aji->second+i_cell]; + U_vector[uji_aji->first + i_cell] = A_vector[uji_aji->second + i_cell]; ++uji_aji; } for (std::size_t i_cell = 0; i_cell < n_cells; ++i_cell) - L_vector[std::get<0>(lii_nuji_nlji)+i_cell] = 1.0; + L_vector[std::get<0>(lii_nuji_nlji) + i_cell] = 1.0; for (std::size_t i = 0; i < std::get<2>(lii_nuji_nlji); ++i) { for (std::size_t i_cell = 0; i_cell < n_cells; ++i_cell) - L_vector[lji_aji->first+i_cell] = A_vector[lji_aji->second+i_cell]; + L_vector[lji_aji->first + i_cell] = A_vector[lji_aji->second + i_cell]; ++lji_aji; } } for (auto& fill_uji : fill_uji_) for (std::size_t i_cell = 0; i_cell < n_cells; ++i_cell) - U_vector[fill_uji+i_cell] = 0; + U_vector[fill_uji + i_cell] = 0; for (auto& fill_lji : fill_lji_) for (std::size_t i_cell = 0; i_cell < n_cells; ++i_cell) - L_vector[fill_lji+i_cell] = 0; + L_vector[fill_lji + i_cell] = 0; for (std::size_t i = 0; i < n; ++i) { for (std::size_t i_cell = 0; i_cell < n_cells; ++i_cell) diff --git a/include/micm/solver/rosenbrock.hpp b/include/micm/solver/rosenbrock.hpp index af7f157be..72dba2f3d 100644 --- a/include/micm/solver/rosenbrock.hpp +++ b/include/micm/solver/rosenbrock.hpp @@ -97,10 +97,10 @@ namespace micm /// @param alpha template void AlphaMinusJacobian(SparseMatrixPolicy& jacobian, const double& alpha) const - requires(!VectorizableSparse); + requires(!VectorizableSparse); template void AlphaMinusJacobian(SparseMatrixPolicy& jacobian, const double& alpha) const - requires(VectorizableSparse); + requires(VectorizableSparse); /// @brief Perform the LU decomposition of the matrix /// @param alpha The alpha value @@ -116,10 +116,10 @@ namespace micm /// @return template double NormalizedError(const DenseMatrixPolicy& y, const DenseMatrixPolicy& y_new, const DenseMatrixPolicy& errors) const - requires(!VectorizableDense); + requires(!VectorizableDense); template double NormalizedError(const DenseMatrixPolicy& y, const DenseMatrixPolicy& y_new, const DenseMatrixPolicy& errors) const - requires(VectorizableDense); + requires(VectorizableDense); }; // end of Abstract Rosenbrock Solver template diff --git a/include/micm/solver/rosenbrock.inl b/include/micm/solver/rosenbrock.inl index dd1b9b9db..cf3af5e70 100644 --- a/include/micm/solver/rosenbrock.inl +++ b/include/micm/solver/rosenbrock.inl @@ -196,7 +196,8 @@ namespace micm template inline void AbstractRosenbrockSolver::AlphaMinusJacobian( SparseMatrixPolicy& jacobian, - const double& alpha) const requires(!VectorizableSparse) + const double& alpha) const + requires(!VectorizableSparse) { MICM_PROFILE_FUNCTION(); @@ -212,7 +213,8 @@ namespace micm template inline void AbstractRosenbrockSolver::AlphaMinusJacobian( SparseMatrixPolicy& jacobian, - const double& alpha) const requires(VectorizableSparse) + const double& alpha) const + requires(VectorizableSparse) { MICM_PROFILE_FUNCTION(); @@ -246,7 +248,8 @@ namespace micm inline double AbstractRosenbrockSolver::NormalizedError( const DenseMatrixPolicy& Y, const DenseMatrixPolicy& Ynew, - const DenseMatrixPolicy& errors) const requires(!VectorizableDense) + const DenseMatrixPolicy& errors) const + requires(!VectorizableDense) { // Solving Ordinary Differential Equations II, page 123 // https://link-springer-com.cuucar.idm.oclc.org/book/10.1007/978-3-642-05221-7 @@ -281,7 +284,8 @@ namespace micm inline double AbstractRosenbrockSolver::NormalizedError( const DenseMatrixPolicy& Y, const DenseMatrixPolicy& Ynew, - const DenseMatrixPolicy& errors) const requires(VectorizableDense) + const DenseMatrixPolicy& errors) const + requires(VectorizableDense) { // Solving Ordinary Differential Equations II, page 123 // https://link-springer-com.cuucar.idm.oclc.org/book/10.1007/978-3-642-05221-7 diff --git a/include/micm/solver/solver.hpp b/include/micm/solver/solver.hpp index f8fc816d0..187d6e46d 100644 --- a/include/micm/solver/solver.hpp +++ b/include/micm/solver/solver.hpp @@ -60,7 +60,9 @@ namespace micm SolverResult Solve(double time_step, StatePolicy& state) { - return solver_.Solve(time_step, state); + auto result = solver_.Solve(time_step, state); + state.variables_.Max(0.0); + return result; } /// @brief Returns the number of grid cells diff --git a/include/micm/solver/solver_result.hpp b/include/micm/solver/solver_result.hpp index fb9797df5..13ea33dc9 100644 --- a/include/micm/solver/solver_result.hpp +++ b/include/micm/solver/solver_result.hpp @@ -48,6 +48,9 @@ namespace micm /// @brief Set all member variables to zero void Reset(); + + /// @brief Print the solver stats to the console + void Print() const; }; inline void SolverStats::Reset() @@ -61,6 +64,17 @@ namespace micm solves_ = 0; } + inline void SolverStats::Print() const + { + std::cout << "Function calls: " << function_calls_ << std::endl; + std::cout << "Jacobian updates: " << jacobian_updates_ << std::endl; + std::cout << "Number of steps: " << number_of_steps_ << std::endl; + std::cout << "Accepted steps: " << accepted_ << std::endl; + std::cout << "Rejected steps: " << rejected_ << std::endl; + std::cout << "Decompositions: " << decompositions_ << std::endl; + std::cout << "Solves: " << solves_ << std::endl; + } + inline std::string SolverStateToString(const SolverState& state) { switch (state) diff --git a/include/micm/solver/state.hpp b/include/micm/solver/state.hpp index 1df4f2783..469faabcf 100644 --- a/include/micm/solver/state.hpp +++ b/include/micm/solver/state.hpp @@ -33,11 +33,12 @@ namespace micm std::set> nonzero_jacobian_elements_{}; }; - template + template< + class DenseMatrixPolicy = StandardDenseMatrix, + class SparseMatrixPolicy = StandardSparseMatrix, + class LuDecompositionPolicy = LuDecomposition, + class LMatrixPolicy = SparseMatrixPolicy, + class UMatrixPolicy = SparseMatrixPolicy> struct State { /// Type of the DenseMatrixPolicy diff --git a/include/micm/solver/state.inl b/include/micm/solver/state.inl index 0d3010484..a1b086c1d 100644 --- a/include/micm/solver/state.inl +++ b/include/micm/solver/state.inl @@ -57,7 +57,12 @@ inline std::error_code make_error_code(MicmStateErrc e) namespace micm { - template + template< + class DenseMatrixPolicy, + class SparseMatrixPolicy, + class LuDecompositionPolicy, + class LMatrixPolicy, + class UMatrixPolicy> inline State::State() : conditions_(), variables_(), @@ -67,8 +72,14 @@ namespace micm { } - template - inline State::State(const StateParameters& parameters) + template< + class DenseMatrixPolicy, + class SparseMatrixPolicy, + class LuDecompositionPolicy, + class LMatrixPolicy, + class UMatrixPolicy> + inline State::State( + const StateParameters& parameters) : conditions_(parameters.number_of_grid_cells_), variables_(parameters.number_of_grid_cells_, parameters.variable_names_.size(), 0.0), custom_rate_parameters_(parameters.number_of_grid_cells_, parameters.custom_rate_parameter_labels_.size(), 0.0), @@ -107,8 +118,14 @@ namespace micm } } - template - inline void State::SetConcentrations( + template< + class DenseMatrixPolicy, + class SparseMatrixPolicy, + class LuDecompositionPolicy, + class LMatrixPolicy, + class UMatrixPolicy> + inline void + State::SetConcentrations( const std::unordered_map>& species_to_concentration) { const std::size_t num_grid_cells = conditions_.size(); @@ -116,8 +133,16 @@ namespace micm SetConcentration({ pair.first }, pair.second); } - template - inline void State::SetConcentration(const Species& species, double concentration) + template< + class DenseMatrixPolicy, + class SparseMatrixPolicy, + class LuDecompositionPolicy, + class LMatrixPolicy, + class UMatrixPolicy> + inline void + State::SetConcentration( + const Species& species, + double concentration) { auto var = variable_map_.find(species.name_); if (var == variable_map_.end()) @@ -127,8 +152,14 @@ namespace micm variables_[0][variable_map_[species.name_]] = concentration; } - template - inline void State::SetConcentration( + template< + class DenseMatrixPolicy, + class SparseMatrixPolicy, + class LuDecompositionPolicy, + class LMatrixPolicy, + class UMatrixPolicy> + inline void + State::SetConcentration( const Species& species, const std::vector& concentration) { @@ -142,9 +173,14 @@ namespace micm variables_[i][i_species] = concentration[i]; } - template - inline void State::UnsafelySetCustomRateParameters( - const std::vector>& parameters) + template< + class DenseMatrixPolicy, + class SparseMatrixPolicy, + class LuDecompositionPolicy, + class LMatrixPolicy, + class UMatrixPolicy> + inline void State:: + UnsafelySetCustomRateParameters(const std::vector>& parameters) { if (parameters.size() != variables_.NumRows()) throw std::system_error( @@ -159,16 +195,30 @@ namespace micm } } - template - inline void State::SetCustomRateParameters( + template< + class DenseMatrixPolicy, + class SparseMatrixPolicy, + class LuDecompositionPolicy, + class LMatrixPolicy, + class UMatrixPolicy> + inline void + State::SetCustomRateParameters( const std::unordered_map>& parameters) { for (auto& pair : parameters) SetCustomRateParameter(pair.first, pair.second); } - template - inline void State::SetCustomRateParameter(const std::string& label, double value) + template< + class DenseMatrixPolicy, + class SparseMatrixPolicy, + class LuDecompositionPolicy, + class LMatrixPolicy, + class UMatrixPolicy> + 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()) @@ -179,8 +229,14 @@ namespace micm custom_rate_parameters_[0][param->second] = value; } - template - inline void State::SetCustomRateParameter( + template< + class DenseMatrixPolicy, + class SparseMatrixPolicy, + class LuDecompositionPolicy, + class LMatrixPolicy, + class UMatrixPolicy> + inline void + State::SetCustomRateParameter( const std::string& label, const std::vector& values) { @@ -194,8 +250,14 @@ namespace micm custom_rate_parameters_[i][param->second] = values[i]; } - template - inline void State::PrintHeader() + template< + class DenseMatrixPolicy, + class SparseMatrixPolicy, + class LuDecompositionPolicy, + class LMatrixPolicy, + class UMatrixPolicy> + inline void + State::PrintHeader() { auto largest_str_iter = std::max_element( variable_names_.begin(), variable_names_.end(), [](const auto& a, const auto& b) { return a.size() < b.size(); }); @@ -215,8 +277,14 @@ namespace micm std::cout << std::endl; } - template - inline void State::PrintState(double time) + template< + class DenseMatrixPolicy, + class SparseMatrixPolicy, + class LuDecompositionPolicy, + class LMatrixPolicy, + class UMatrixPolicy> + inline void State::PrintState( + double time) { std::ios oldState(nullptr); oldState.copyfmt(std::cout); diff --git a/include/micm/util/constants.hpp b/include/micm/util/constants.hpp index 20f607466..0a4217366 100644 --- a/include/micm/util/constants.hpp +++ b/include/micm/util/constants.hpp @@ -9,5 +9,5 @@ namespace micm static constexpr double BOLTZMANN_CONSTANT = 1.380649e-23; // J K^{-1} static constexpr double AVOGADRO_CONSTANT = 6.02214076e23; // # mol^{-1} static constexpr double GAS_CONSTANT = BOLTZMANN_CONSTANT * AVOGADRO_CONSTANT; // J K^{-1} mol^{-1} - } // namespace constants + } // namespace constants } // namespace micm diff --git a/include/micm/util/matrix.hpp b/include/micm/util/matrix.hpp index 5c55e50ec..bd63456ac 100644 --- a/include/micm/util/matrix.hpp +++ b/include/micm/util/matrix.hpp @@ -15,8 +15,7 @@ namespace micm /// Concept for vectorizable matrices template - concept VectorizableDense = requires(T t) - { + concept VectorizableDense = requires(T t) { t.GroupSize(); t.GroupVectorSize(); t.NumberOfGroups(); @@ -231,6 +230,22 @@ namespace micm y += alpha * (*(x_iter++)); } + /// @brief For each element of the matrix, perform y = max(y, x), where x is a scalar constant + /// @param x The scalar constant to compare against + void Max(const T &x) + { + for (auto &y : data_) + y = std::max(y, x); + } + + /// @brief For each element of the matrix, perform y = min(y, x), where x is a scalar constant + /// @param x The scalar constant to compare against + void Min(const T &x) + { + for (auto &y : data_) + y = std::min(y, x); + } + void ForEach(const std::function f, const Matrix &a) { auto a_iter = a.AsVector().begin(); diff --git a/include/micm/util/sparse_matrix.hpp b/include/micm/util/sparse_matrix.hpp index 81e42e3fe..bae5239b4 100644 --- a/include/micm/util/sparse_matrix.hpp +++ b/include/micm/util/sparse_matrix.hpp @@ -17,16 +17,14 @@ namespace micm { /// Concept for vectorizable matrices template - concept VectorizableSparse = requires(T t) - { + concept VectorizableSparse = requires(T t) { t.GroupSize(); t.GroupVectorSize(); t.NumberOfGroups(0); }; template - concept SparseMatrixConcept = requires(T t) - { + concept SparseMatrixConcept = requires(T t) { t.NumRows(); t.NumColumns(); t.NumberOfBlocks(); @@ -55,8 +53,8 @@ namespace micm using value_type = T; protected: - std::size_t number_of_blocks_; // Number of block sub-matrices in the overall matrix - std::size_t block_size_; // Size of each block sub-matrix (number of rows or columns per block) + std::size_t number_of_blocks_; // Number of block sub-matrices in the overall matrix + std::size_t block_size_; // Size of each block sub-matrix (number of rows or columns per block) std::size_t number_of_non_zero_elements_per_block_; // Number of non-zero elements in each block std::vector data_; // Value of each non-zero matrix element diff --git a/include/micm/util/sparse_matrix_standard_ordering_compressed_sparse_column.hpp b/include/micm/util/sparse_matrix_standard_ordering_compressed_sparse_column.hpp index 7ab383c49..1dedc1179 100644 --- a/include/micm/util/sparse_matrix_standard_ordering_compressed_sparse_column.hpp +++ b/include/micm/util/sparse_matrix_standard_ordering_compressed_sparse_column.hpp @@ -6,12 +6,12 @@ #include #include -#include #include -#include -#include #include +#include #include +#include +#include namespace micm { @@ -65,11 +65,7 @@ namespace micm /// @param row Index of the row in the block /// @param column Index of the column in the block /// @return Index of the element in the compressed data vector - std::size_t VectorIndex( - std::size_t number_of_blocks, - std::size_t block, - std::size_t row, - std::size_t column) const + std::size_t VectorIndex(std::size_t number_of_blocks, std::size_t block, std::size_t row, std::size_t column) const { if (column >= column_start_.size() - 1 || row >= column_start_.size() - 1 || block >= number_of_blocks) throw std::system_error(make_error_code(MicmMatrixErrc::ElementOutOfRange)); @@ -85,12 +81,10 @@ namespace micm /// @param number_of_blocks Total number of block sub-matrices in the overall matrix /// @param data Compressed data vector /// @param value Value to add to the diagonal - void AddToDiagonal( - const std::size_t number_of_blocks, - auto& data, - auto value) const + void AddToDiagonal(const std::size_t number_of_blocks, auto& data, auto value) const { - for (std::size_t block_start = 0; block_start < number_of_blocks * column_ids_.size(); block_start += column_ids_.size()) + for (std::size_t block_start = 0; block_start < number_of_blocks * column_ids_.size(); + block_start += column_ids_.size()) for (const auto& i : diagonal_ids_) data[block_start + i] += value; } @@ -110,11 +104,12 @@ namespace micm } private: - /// @brief Returns the column ids of each non-zero element in a block /// @param block_size Number of rows or columns for each block /// @param non_zero_elements Set of non-zero elements in the matrix - static std::vector ColumnIdsVector(const std::size_t block_size, const std::set> non_zero_elements) + static std::vector ColumnIdsVector( + const std::size_t block_size, + const std::set> non_zero_elements) { std::vector ids; ids.reserve(non_zero_elements.size()); @@ -132,7 +127,9 @@ namespace micm /// @brief Returns the start and end indices of each column in a block in column_ids_ /// @param block_size Number of rows or columns for each block /// @param non_zero_elements Set of non-zero elements in the matrix - static std::vector ColumnStartVector(const std::size_t block_size, const std::set> non_zero_elements) + static std::vector ColumnStartVector( + const std::size_t block_size, + const std::set> non_zero_elements) { std::vector starts(block_size + 1, 0); std::size_t total_elem = 0; @@ -151,7 +148,6 @@ namespace micm } public: - /// @brief Returns whether a particular element is always zero /// @param row Row index /// @param column Column index @@ -166,4 +162,4 @@ namespace micm return (elem == end); } }; -} // namespace micm \ No newline at end of file +} // namespace micm \ No newline at end of file diff --git a/include/micm/util/sparse_matrix_standard_ordering_compressed_sparse_row.hpp b/include/micm/util/sparse_matrix_standard_ordering_compressed_sparse_row.hpp index 5b17e49f2..60a28b98a 100644 --- a/include/micm/util/sparse_matrix_standard_ordering_compressed_sparse_row.hpp +++ b/include/micm/util/sparse_matrix_standard_ordering_compressed_sparse_row.hpp @@ -6,12 +6,12 @@ #include #include -#include #include -#include -#include #include +#include #include +#include +#include namespace micm { @@ -23,7 +23,7 @@ namespace micm class SparseMatrixStandardOrderingCompressedSparseRow { /// @brief Row ids of each non-zero element in a block - std::vector row_ids_; + std::vector row_ids_; /// @brief Start and end indices of each row in a block in row_ids_ std::vector row_start_; /// @brief Indices along the diagonal of each block @@ -65,11 +65,7 @@ namespace micm /// @param row Index of the row in the block /// @param column Index of the column in the block /// @return Index of the element in the compressed data vector - std::size_t VectorIndex( - std::size_t number_of_blocks, - std::size_t block, - std::size_t row, - std::size_t column) const + std::size_t VectorIndex(std::size_t number_of_blocks, std::size_t block, std::size_t row, std::size_t column) const { if (row >= row_start_.size() - 1 || column >= row_start_.size() - 1 || block >= number_of_blocks) throw std::system_error(make_error_code(MicmMatrixErrc::ElementOutOfRange)); @@ -85,10 +81,7 @@ namespace micm /// @param number_of_blocks Total number of block sub-matrices in the overall matrix /// @param data Compressed data vector /// @param value Value to add to the diagonal - void AddToDiagonal( - const std::size_t number_of_blocks, - auto& data, - auto value) const + void AddToDiagonal(const std::size_t number_of_blocks, auto& data, auto value) const { for (std::size_t block_start = 0; block_start < number_of_blocks * row_ids_.size(); block_start += row_ids_.size()) for (const auto& i : diagonal_ids_) @@ -109,12 +102,13 @@ namespace micm return indices; } - private: - + private: /// @brief Returns the row ids of each non-zero element in a block /// @param block_size Number of rows or columns for each block /// @param non_zero_elements Set of non-zero elements in the matrix - static std::vector RowIdsVector(const std::size_t block_size, const std::set> non_zero_elements) + static std::vector RowIdsVector( + const std::size_t block_size, + const std::set> non_zero_elements) { std::vector ids; ids.reserve(non_zero_elements.size()); @@ -129,7 +123,9 @@ namespace micm /// @brief Returns the start and end indices of each row in a block in row_ids_ /// @param block_size Number of rows or columns for each block /// @param non_zero_elements Set of non-zero elements in the matrix - static std::vector RowStartVector(const std::size_t block_size, const std::set> non_zero_elements) + static std::vector RowStartVector( + const std::size_t block_size, + const std::set> non_zero_elements) { std::vector starts(block_size + 1, 0); std::size_t total_elem = 0; @@ -144,8 +140,7 @@ namespace micm return starts; } - public: - + public: /// @brief Returns whether a particular element is always zero /// @param row Row index /// @param column Column index diff --git a/include/micm/util/sparse_matrix_vector_ordering.hpp b/include/micm/util/sparse_matrix_vector_ordering.hpp index 4daa0bd09..6bdf6cff0 100644 --- a/include/micm/util/sparse_matrix_vector_ordering.hpp +++ b/include/micm/util/sparse_matrix_vector_ordering.hpp @@ -2,8 +2,8 @@ // SPDX-License-Identifier: Apache-2.0 #pragma once -#include "sparse_matrix_vector_ordering_compressed_sparse_row.hpp" #include "sparse_matrix_vector_ordering_compressed_sparse_column.hpp" +#include "sparse_matrix_vector_ordering_compressed_sparse_row.hpp" namespace micm { diff --git a/include/micm/util/sparse_matrix_vector_ordering_compressed_sparse_column.hpp b/include/micm/util/sparse_matrix_vector_ordering_compressed_sparse_column.hpp index 3a21fc9e8..a77cac4b6 100644 --- a/include/micm/util/sparse_matrix_vector_ordering_compressed_sparse_column.hpp +++ b/include/micm/util/sparse_matrix_vector_ordering_compressed_sparse_column.hpp @@ -11,9 +11,9 @@ #include #include #include +#include #include #include -#include namespace micm { @@ -64,7 +64,7 @@ namespace micm /// @return Size of the compressed data vector std::size_t VectorSize(std::size_t number_of_blocks) const { - return std::ceil((double)number_of_blocks / (double)L) * L *column_ids_.size(); + return std::ceil((double)number_of_blocks / (double)L) * L * column_ids_.size(); } /// @brief Returns the index for a particular element in the compressed data vector @@ -73,11 +73,7 @@ namespace micm /// @param row Index of the row in the block /// @param column Index of the column in the block /// @return Index of the element in the compressed data vector - std::size_t VectorIndex( - std::size_t number_of_blocks, - std::size_t block, - std::size_t row, - std::size_t column) const + std::size_t VectorIndex(std::size_t number_of_blocks, std::size_t block, std::size_t row, std::size_t column) const { if (column >= column_start_.size() - 1 || row >= column_start_.size() - 1 || block >= number_of_blocks) throw std::system_error(make_error_code(MicmMatrixErrc::ElementOutOfRange)); @@ -93,10 +89,7 @@ namespace micm /// @param number_of_blocks Total number of block sub-matrices in the overall matrix /// @param data Compressed data vector /// @param value Value to add to the diagonal - void AddToDiagonal( - const std::size_t number_of_blocks, - auto& data, - const auto value) const + void AddToDiagonal(const std::size_t number_of_blocks, auto& data, const auto value) const { for (std::size_t i_group = 0; i_group < number_of_blocks; i_group += L) { @@ -149,7 +142,6 @@ namespace micm } private: - /// @brief Returns the column ids of each non-zero element in a block /// @param block_size Number of rows or columns in each block /// @param non_zero_elements Set of non-zero elements in the matrix @@ -192,11 +184,10 @@ namespace micm ++total_elem; } starts[curr_row + 1] = total_elem; - return starts; + return starts; } public: - /// @brief Returns whether a given row and column index is a zero element /// @param row Index of the row /// @param column Index of the column @@ -211,4 +202,4 @@ namespace micm return (elem == end); } }; -} \ No newline at end of file +} // namespace micm \ No newline at end of file diff --git a/include/micm/util/sparse_matrix_vector_ordering_compressed_sparse_row.hpp b/include/micm/util/sparse_matrix_vector_ordering_compressed_sparse_row.hpp index 609b3bea8..08bac1800 100644 --- a/include/micm/util/sparse_matrix_vector_ordering_compressed_sparse_row.hpp +++ b/include/micm/util/sparse_matrix_vector_ordering_compressed_sparse_row.hpp @@ -11,9 +11,9 @@ #include #include #include +#include #include #include -#include namespace micm { @@ -73,11 +73,7 @@ namespace micm /// @param row Index of the row in the block /// @param column Index of the column in the block /// @return Index of the element in the compressed data vector - std::size_t VectorIndex( - std::size_t number_of_blocks, - std::size_t block, - std::size_t row, - std::size_t column) const + std::size_t VectorIndex(std::size_t number_of_blocks, std::size_t block, std::size_t row, std::size_t column) const { if (row >= row_start_.size() - 1 || column >= row_start_.size() - 1 || block >= number_of_blocks) throw std::system_error(make_error_code(MicmMatrixErrc::ElementOutOfRange)); @@ -93,10 +89,7 @@ namespace micm /// @param number_of_blocks Total number of block sub-matrices in the overall matrix /// @param data Compressed data vector /// @param value Value to add to the diagonal - void AddToDiagonal( - const std::size_t number_of_blocks, - auto& data, - auto value) const + void AddToDiagonal(const std::size_t number_of_blocks, auto& data, auto value) const { for (std::size_t i_group = 0; i_group < number_of_blocks; i_group += L) { @@ -149,11 +142,12 @@ namespace micm } private: - /// @brief Returns the row ids of each non-zero element in a block /// @param block_size Number of rows or columns for each block /// @param non_zero_elements Set of non-zero elements in the matrix - static std::vector RowIdsVector(const std::size_t block_size, const std::set> non_zero_elements) + static std::vector RowIdsVector( + const std::size_t block_size, + const std::set> non_zero_elements) { std::vector ids; ids.reserve(non_zero_elements.size()); @@ -168,7 +162,9 @@ namespace micm /// @brief Returns the start and end indices of each row in a block in row_ids_ /// @param block_size Number of rows or columns for each block /// @param non_zero_elements Set of non-zero elements in the matrix - static std::vector RowStartVector(const std::size_t block_size, const std::set> non_zero_elements) + static std::vector RowStartVector( + const std::size_t block_size, + const std::set> non_zero_elements) { std::vector starts(block_size + 1, 0); std::size_t total_elem = 0; @@ -183,8 +179,7 @@ namespace micm return starts; } - public: - + public: /// @brief Returns whether a particular element is always zero /// @param row Row index /// @param column Column index diff --git a/include/micm/util/vector_matrix.hpp b/include/micm/util/vector_matrix.hpp index 0ef022f30..ecc63cb47 100644 --- a/include/micm/util/vector_matrix.hpp +++ b/include/micm/util/vector_matrix.hpp @@ -288,6 +288,26 @@ namespace micm } } + /// @brief For each element of the VectorMatrix, perform y = max(y, x), where x is a scalar constant + /// @param x The scalar constant to compare against + void Max(const T &x) + { + MICM_PROFILE_FUNCTION(); + + for (auto &y : data_) + y = std::max(y, x); + } + + /// @brief For each element of the VectorMatrix, perform y = min(y, x), where x is a scalar constant + /// @param x The scalar constant to compare against + void Min(const T &x) + { + MICM_PROFILE_FUNCTION(); + + for (auto &y : data_) + y = std::min(y, x); + } + void ForEach(const std::function f, const VectorMatrix &a) { auto this_iter = data_.begin(); diff --git a/src/process/process_set.cu b/src/process/process_set.cu index 9f36a4afa..78ed5bd65 100644 --- a/src/process/process_set.cu +++ b/src/process/process_set.cu @@ -56,8 +56,8 @@ namespace micm prod_id_offset += d_number_of_products[i_rxn]; yield_offset += d_number_of_products[i_rxn]; } // end of loop over number of reactions - } // end of checking a valid CUDA thread id - } // end of AddForcingTerms_kernel + } // end of checking a valid CUDA thread id + } // end of AddForcingTerms_kernel /// This is the CUDA kernel that forms the negative Jacobian matrix (-J) on the device __global__ void SubtractJacobianTermsKernel( @@ -117,8 +117,8 @@ namespace micm react_ids_offset += d_number_of_reactants[i_rxn]; yields_offset += d_number_of_products[i_rxn]; } // end of loop over reactions in a grid cell - } // end of checking a CUDA thread id - } // end of SubtractJacobianTermsKernel + } // end of checking a CUDA thread id + } // end of SubtractJacobianTermsKernel /// This is the function that will copy the constant data /// members of class "CudaProcessSet" to the device, @@ -303,5 +303,5 @@ namespace micm micm::cuda::CudaStreamSingleton::GetInstance().GetCudaStream(0)>>>( rate_constants_param, state_variables_param, forcing_param, devstruct); } // end of AddForcingTermsKernelDriver - } // namespace cuda + } // namespace cuda } // namespace micm diff --git a/src/solver/lu_decomposition_doolittle.cu b/src/solver/lu_decomposition_doolittle.cu index e942a16fd..cc826589a 100644 --- a/src/solver/lu_decomposition_doolittle.cu +++ b/src/solver/lu_decomposition_doolittle.cu @@ -295,5 +295,5 @@ namespace micm DecomposeKernel<<>>( A_param, L_param, U_param, devstruct); } // end of DecomposeKernelDriver - } // end of namespace cuda + } // end of namespace cuda } // end of namespace micm diff --git a/src/solver/rosenbrock.cu b/src/solver/rosenbrock.cu index 06d355139..e2a980e82 100644 --- a/src/solver/rosenbrock.cu +++ b/src/solver/rosenbrock.cu @@ -380,5 +380,5 @@ namespace micm } // end of if-else for CUDA/CUBLAS implementation return std::max(normalized_error, 1.0e-10); } // end of NormalizedErrorDriver function - } // namespace cuda + } // namespace cuda } // namespace micm diff --git a/src/util/cuda_matrix.cu b/src/util/cuda_matrix.cu index 08a3946bf..97fa1624a 100644 --- a/src/util/cuda_matrix.cu +++ b/src/util/cuda_matrix.cu @@ -35,6 +35,46 @@ namespace micm return err; } + template + __global__ void MatrixMaxKernel(T* d_data, std::size_t number_of_elements, T val) + { + std::size_t tid = blockIdx.x * BLOCK_SIZE + threadIdx.x; + if (tid < number_of_elements) + { + d_data[tid] = d_data[tid] > val ? d_data[tid] : val; + } + } + + template + cudaError_t MatrixMax(CudaMatrixParam& param, T val) + { + std::size_t number_of_blocks = (param.number_of_elements_ + BLOCK_SIZE - 1) / BLOCK_SIZE; + MatrixMaxKernel<<>>( + param.d_data_, param.number_of_elements_, val); + cudaError_t err = cudaGetLastError(); + return err; + } + + template + __global__ void MatrixMinKernel(T* d_data, std::size_t number_of_elements, T val) + { + std::size_t tid = blockIdx.x * BLOCK_SIZE + threadIdx.x; + if (tid < number_of_elements) + { + d_data[tid] = d_data[tid] < val ? d_data[tid] : val; + } + } + + template + cudaError_t MatrixMin(CudaMatrixParam& param, T val) + { + std::size_t number_of_blocks = (param.number_of_elements_ + BLOCK_SIZE - 1) / BLOCK_SIZE; + MatrixMinKernel<<>>( + param.d_data_, param.number_of_elements_, val); + cudaError_t err = cudaGetLastError(); + return err; + } + template cudaError_t CopyToDevice(CudaMatrixParam& param, std::vector& h_data) { @@ -98,6 +138,8 @@ namespace micm // source code needs the instantiation of the template template cudaError_t MallocVector(CudaMatrixParam& param, std::size_t number_of_elements); template cudaError_t MallocVector(CudaMatrixParam& param, std::size_t number_of_elements); + template cudaError_t MatrixMax(CudaMatrixParam& param, double val); + template cudaError_t MatrixMin(CudaMatrixParam& param, double val); template cudaError_t CopyToDevice(CudaMatrixParam& param, std::vector& h_data); template cudaError_t CopyToHost(CudaMatrixParam& param, std::vector& h_data); template cudaError_t CopyToDeviceFromDevice( diff --git a/test/integration/test_analytical_backward_euler.cpp b/test/integration/test_analytical_backward_euler.cpp index 7a17e56f5..b4a0c2850 100644 --- a/test/integration/test_analytical_backward_euler.cpp +++ b/test/integration/test_analytical_backward_euler.cpp @@ -39,8 +39,10 @@ using VectorBackwardEulerDolittleCSC = micm::CpuSolverBuilder< micm::LuDecompositionDoolittle>; template -using VectorStateTypeDoolittleCSC = - micm::State, micm::SparseMatrix>, micm::LuDecompositionDoolittle>; +using VectorStateTypeDoolittleCSC = micm::State< + micm::VectorMatrix, + micm::SparseMatrix>, + micm::LuDecompositionDoolittle>; template using VectorBackwardEulerMozart = micm::CpuSolverBuilder< @@ -59,11 +61,33 @@ template using VectorBackwardEulerMozartCSC = micm::CpuSolverBuilder< micm::BackwardEulerSolverParameters, micm::VectorMatrix, - micm::SparseMatrix>, micm::LuDecompositionMozart>; + micm::SparseMatrix>, + micm::LuDecompositionMozart>; + +template +using VectorStateTypeMozartCSC = micm::State< + micm::VectorMatrix, + micm::SparseMatrix>, + micm::LuDecompositionMozart>; template -using VectorStateTypeMozartCSC = - micm::State, micm::SparseMatrix>, micm::LuDecompositionMozart>; +using VectorBackwardEulerMozartInPlace = micm::SolverBuilder< + micm::BackwardEulerSolverParameters, + micm::VectorMatrix, + micm::SparseMatrix>, + micm::ProcessSet, + micm::LuDecompositionMozartInPlace, + micm::LinearSolverInPlace>, micm::LuDecompositionMozartInPlace>, + micm::State< + micm::VectorMatrix, + micm::SparseMatrix>, + micm::LuDecompositionMozartInPlace>>; + +template +using VectorStateTypeMozartInPlace = micm::State< + micm::VectorMatrix, + micm::SparseMatrix>, + micm::LuDecompositionMozartInPlace>; template using VectorBackwardEulerMozartInPlace = micm::SolverBuilder< @@ -126,18 +150,26 @@ TEST(AnalyticalExamples, Troe) backward_euler_vector_doolittle_3, 1e-6); test_analytical_troe, VectorStateTypeDoolittle<4>>( backward_euler_vector_doolittle_4, 1e-6); - test_analytical_troe, VectorStateTypeDoolittleCSC<1>>(backward_euler_vector_doolittle_csc_1, 1e-6); - test_analytical_troe, VectorStateTypeDoolittleCSC<2>>(backward_euler_vector_doolittle_csc_2, 1e-6); - test_analytical_troe, VectorStateTypeDoolittleCSC<3>>(backward_euler_vector_doolittle_csc_3, 1e-6); - test_analytical_troe, VectorStateTypeDoolittleCSC<4>>(backward_euler_vector_doolittle_csc_4, 1e-6); + test_analytical_troe, VectorStateTypeDoolittleCSC<1>>( + backward_euler_vector_doolittle_csc_1, 1e-6); + test_analytical_troe, VectorStateTypeDoolittleCSC<2>>( + backward_euler_vector_doolittle_csc_2, 1e-6); + test_analytical_troe, VectorStateTypeDoolittleCSC<3>>( + backward_euler_vector_doolittle_csc_3, 1e-6); + test_analytical_troe, VectorStateTypeDoolittleCSC<4>>( + backward_euler_vector_doolittle_csc_4, 1e-6); test_analytical_troe, VectorStateTypeMozart<1>>(backward_euler_vector_mozart_1, 1e-6); test_analytical_troe, VectorStateTypeMozart<2>>(backward_euler_vector_mozart_2, 1e-6); test_analytical_troe, VectorStateTypeMozart<3>>(backward_euler_vector_mozart_3, 1e-6); test_analytical_troe, VectorStateTypeMozart<4>>(backward_euler_vector_mozart_4, 1e-6); - test_analytical_troe, VectorStateTypeMozartCSC<1>>(backward_euler_vector_mozart_csc_1, 1e-6); - test_analytical_troe, VectorStateTypeMozartCSC<2>>(backward_euler_vector_mozart_csc_2, 1e-6); - test_analytical_troe, VectorStateTypeMozartCSC<3>>(backward_euler_vector_mozart_csc_3, 1e-6); - test_analytical_troe, VectorStateTypeMozartCSC<4>>(backward_euler_vector_mozart_csc_4, 1e-6); + test_analytical_troe, VectorStateTypeMozartCSC<1>>( + backward_euler_vector_mozart_csc_1, 1e-6); + test_analytical_troe, VectorStateTypeMozartCSC<2>>( + backward_euler_vector_mozart_csc_2, 1e-6); + test_analytical_troe, VectorStateTypeMozartCSC<3>>( + backward_euler_vector_mozart_csc_3, 1e-6); + test_analytical_troe, VectorStateTypeMozartCSC<4>>( + backward_euler_vector_mozart_csc_4, 1e-6); test_analytical_troe, VectorStateTypeMozartInPlace<1>>( backward_euler_vector_mozart_in_place_1, 1e-6); test_analytical_troe, VectorStateTypeMozartInPlace<2>>( diff --git a/test/tutorial/test_multiple_grid_cells.cpp b/test/tutorial/test_multiple_grid_cells.cpp index 173a62440..d011919ba 100644 --- a/test/tutorial/test_multiple_grid_cells.cpp +++ b/test/tutorial/test_multiple_grid_cells.cpp @@ -86,7 +86,7 @@ int main() { solver.CalculateRateConstants(state); auto result = solver.Solve(time_step - elapsed_solve_time, state); - elapsed_solve_time = result.final_time_; + elapsed_solve_time += result.final_time_; } state.PrintState(time_step * (i + 1)); } diff --git a/test/tutorial/test_rate_constants_no_user_defined_by_hand.cpp b/test/tutorial/test_rate_constants_no_user_defined_by_hand.cpp index f105da91b..7172e13f1 100644 --- a/test/tutorial/test_rate_constants_no_user_defined_by_hand.cpp +++ b/test/tutorial/test_rate_constants_no_user_defined_by_hand.cpp @@ -138,7 +138,7 @@ int main(const int argc, const char* argv[]) while (elapsed_solve_time < time_step) { auto result = solver.Solve(time_step - elapsed_solve_time, state); - elapsed_solve_time = result.final_time_; + elapsed_solve_time += result.final_time_; } state.PrintState(time_step * (i + 1)); diff --git a/test/tutorial/test_rate_constants_no_user_defined_with_config.cpp b/test/tutorial/test_rate_constants_no_user_defined_with_config.cpp index 873c7f28b..1b65cfd3a 100644 --- a/test/tutorial/test_rate_constants_no_user_defined_with_config.cpp +++ b/test/tutorial/test_rate_constants_no_user_defined_with_config.cpp @@ -80,7 +80,7 @@ int main(const int argc, const char* argv[]) while (elapsed_solve_time < time_step) { auto result = solver.Solve(time_step - elapsed_solve_time, state); - elapsed_solve_time = result.final_time_; + elapsed_solve_time += result.final_time_; } state.PrintState(time_step * (i + 1)); diff --git a/test/tutorial/test_rate_constants_user_defined_by_hand.cpp b/test/tutorial/test_rate_constants_user_defined_by_hand.cpp index 7ded8f81a..8475edecc 100644 --- a/test/tutorial/test_rate_constants_user_defined_by_hand.cpp +++ b/test/tutorial/test_rate_constants_user_defined_by_hand.cpp @@ -164,7 +164,7 @@ int main(const int argc, const char* argv[]) while (elapsed_solve_time < time_step) { auto result = solver.Solve(time_step - elapsed_solve_time, state); - elapsed_solve_time = result.final_time_; + elapsed_solve_time += result.final_time_; } state.PrintState(time_step * (i + 1)); diff --git a/test/tutorial/test_rate_constants_user_defined_with_config.cpp b/test/tutorial/test_rate_constants_user_defined_with_config.cpp index 7421d645c..5603a699a 100644 --- a/test/tutorial/test_rate_constants_user_defined_with_config.cpp +++ b/test/tutorial/test_rate_constants_user_defined_with_config.cpp @@ -89,7 +89,7 @@ int main(const int argc, const char* argv[]) while (elapsed_solve_time < time_step) { auto result = solver.Solve(time_step - elapsed_solve_time, state); - elapsed_solve_time = result.final_time_; + elapsed_solve_time += result.final_time_; } state.PrintState(time_step * (i + 1)); diff --git a/test/tutorial/test_solver_configuration.cpp b/test/tutorial/test_solver_configuration.cpp index 522908d23..06ad83129 100644 --- a/test/tutorial/test_solver_configuration.cpp +++ b/test/tutorial/test_solver_configuration.cpp @@ -65,7 +65,7 @@ void test_solver_type(auto& solver) total_stats.decompositions_ += result.stats_.decompositions_; total_stats.solves_ += result.stats_.solves_; - elapsed_solve_time = result.final_time_; + elapsed_solve_time += result.final_time_; } state.PrintState(time_step * (i + 1)); diff --git a/test/tutorial/test_vectorized_matrix_solver.cpp b/test/tutorial/test_vectorized_matrix_solver.cpp index 29538e7d5..62ce9a5d6 100644 --- a/test/tutorial/test_vectorized_matrix_solver.cpp +++ b/test/tutorial/test_vectorized_matrix_solver.cpp @@ -43,7 +43,7 @@ void solve(auto& solver, auto& state, std::size_t number_of_grid_cells) while (elapsed_solve_time < time_step && solver_state != SolverState::Converged) { auto result = solver.Solve(time_step - elapsed_solve_time, state); - elapsed_solve_time = result.final_time_; + elapsed_solve_time += result.final_time_; solver_state = result.state_; } } diff --git a/test/unit/cuda/util/test_cuda_dense_matrix.cpp b/test/unit/cuda/util/test_cuda_dense_matrix.cpp index 8bac62627..0c2f7215d 100644 --- a/test/unit/cuda/util/test_cuda_dense_matrix.cpp +++ b/test/unit/cuda/util/test_cuda_dense_matrix.cpp @@ -571,4 +571,62 @@ TEST(CudaDenseMatrix, CopyFunction) EXPECT_EQ(matrix2[i][j], matrix[i][j]); } } +} + +TEST(CudaDenseMatrix, TestMax) +{ + micm::CudaDenseMatrix matrix{ 2, 3, 0.0 }; + matrix.CopyToDevice(); + matrix.Max(2.0); + matrix.CopyToHost(); + + for (auto& elem : matrix.AsVector()) + { + EXPECT_EQ(elem, 2.0); + } + + for (auto& elem : matrix.AsVector()) + { + elem = 1.0; + } + matrix[1][1] = 3.0; + matrix.CopyToDevice(); + matrix.Max(2.0); + matrix.CopyToHost(); + + EXPECT_EQ(matrix[0][0], 2.0); + EXPECT_EQ(matrix[0][1], 2.0); + EXPECT_EQ(matrix[0][2], 2.0); + EXPECT_EQ(matrix[1][0], 2.0); + EXPECT_EQ(matrix[1][1], 3.0); + EXPECT_EQ(matrix[1][2], 2.0); +} + +TEST(CudaDenseMatrix, TestMin) +{ + micm::CudaDenseMatrix matrix{ 2, 3, 0.0 }; + matrix.CopyToDevice(); + matrix.Min(2.0); + matrix.CopyToHost(); + + for (auto& elem : matrix.AsVector()) + { + EXPECT_EQ(elem, 0.0); + } + + for (auto& elem : matrix.AsVector()) + { + elem = 1.0; + } + matrix[1][1] = 3.0; + matrix.CopyToDevice(); + matrix.Min(2.0); + matrix.CopyToHost(); + + EXPECT_EQ(matrix[0][0], 1.0); + EXPECT_EQ(matrix[0][1], 1.0); + EXPECT_EQ(matrix[0][2], 1.0); + EXPECT_EQ(matrix[1][0], 1.0); + EXPECT_EQ(matrix[1][1], 2.0); + EXPECT_EQ(matrix[1][2], 1.0); } \ No newline at end of file diff --git a/test/unit/solver/test_linear_solver_policy.hpp b/test/unit/solver/test_linear_solver_policy.hpp index 2a632050a..3b6e5d144 100644 --- a/test/unit/solver/test_linear_solver_policy.hpp +++ b/test/unit/solver/test_linear_solver_policy.hpp @@ -11,9 +11,7 @@ template void CopyToDeviceDense(MatrixPolicy& matrix) { if constexpr (requires { - { - matrix.CopyToDevice() - } -> std::same_as; + { matrix.CopyToDevice() } -> std::same_as; }) matrix.CopyToDevice(); } @@ -22,9 +20,7 @@ template void CopyToDeviceSparse(SparseMatrixPolicy& matrix) { if constexpr (requires { - { - matrix.CopyToDevice() - } -> std::same_as; + { matrix.CopyToDevice() } -> std::same_as; }) matrix.CopyToDevice(); } @@ -33,9 +29,7 @@ template void CopyToHostDense(MatrixPolicy& matrix) { if constexpr (requires { - { - matrix.CopyToHost() - } -> std::same_as; + { matrix.CopyToHost() } -> std::same_as; }) matrix.CopyToHost(); } @@ -75,8 +69,7 @@ void print_matrix(const SparseMatrixPolicy& matrix, std::size_t width) { if (matrix.IsZero(i, j)) { - std::cout << " " << std::setfill('-') << std::setw(width) << "-" - << " "; + std::cout << " " << std::setfill('-') << std::setw(width) << "-" << " "; } else { @@ -248,7 +241,8 @@ void testExtremeInitialValue(std::size_t number_of_blocks, double initial_value) CopyToDeviceDense(x); LinearSolverPolicy solver = LinearSolverPolicy(A, initial_value); - auto lu = micm::LuDecomposition::GetLUMatrices(A, initial_value); + auto lu = + micm::LuDecomposition::GetLUMatrices(A, initial_value); auto lower_matrix = std::move(lu.first); auto upper_matrix = std::move(lu.second); @@ -347,10 +341,13 @@ void testMarkowitzReordering() SparseMatrixPolicy reordered_jac{ builder }; auto orig_LU_calc = micm::LuDecomposition::Create(orig_jac); - auto reordered_LU_calc = micm::LuDecomposition::Create(reordered_jac); + auto reordered_LU_calc = + micm::LuDecomposition::Create(reordered_jac); - auto orig_LU = orig_LU_calc.template GetLUMatrices(orig_jac, 0.0); - auto reordered_LU = reordered_LU_calc.template GetLUMatrices(reordered_jac, 0.0); + auto orig_LU = + orig_LU_calc.template GetLUMatrices(orig_jac, 0.0); + auto reordered_LU = reordered_LU_calc.template GetLUMatrices( + reordered_jac, 0.0); std::size_t sum_orig = 0; std::size_t sum_reordered = 0; diff --git a/test/unit/solver/test_lu_decomposition_policy.hpp b/test/unit/solver/test_lu_decomposition_policy.hpp index 2f9001eb2..afea2f18f 100644 --- a/test/unit/solver/test_lu_decomposition_policy.hpp +++ b/test/unit/solver/test_lu_decomposition_policy.hpp @@ -10,9 +10,7 @@ template void CopyToDevice(MatrixPolicy& matrix) { if constexpr (requires { - { - matrix.CopyToDevice() - } -> std::same_as; + { matrix.CopyToDevice() } -> std::same_as; }) matrix.CopyToDevice(); } @@ -21,9 +19,7 @@ template void CopyToHost(MatrixPolicy& matrix) { if constexpr (requires { - { - matrix.CopyToHost() - } -> std::same_as; + { matrix.CopyToHost() } -> std::same_as; }) matrix.CopyToHost(); } @@ -78,8 +74,7 @@ void print_matrix(const auto& matrix, std::size_t width) { if (matrix.IsZero(i, j)) { - std::cout << " " << std::setfill('-') << std::setw(width) << "-" - << " "; + std::cout << " " << std::setfill('-') << std::setw(width) << "-" << " "; } else { @@ -118,7 +113,8 @@ void testDenseMatrix() A[0][2][1] = -2; A[0][2][2] = 8; - LuDecompositionPolicy lud = LuDecompositionPolicy::template Create(A); + LuDecompositionPolicy lud = + LuDecompositionPolicy::template Create(A); auto LU = LuDecompositionPolicy::template GetLUMatrices(A, 0); lud.template Decompose(A, LU.first, LU.second); check_results( @@ -145,7 +141,8 @@ void testRandomMatrix(std::size_t number_of_blocks) for (std::size_t i_block = 0; i_block < number_of_blocks; ++i_block) A[i_block][i][j] = get_double(); - LuDecompositionPolicy lud = LuDecompositionPolicy::template Create(A); + LuDecompositionPolicy lud = + LuDecompositionPolicy::template Create(A); auto LU = LuDecompositionPolicy::template GetLUMatrices(A, 0); lud.template Decompose(A, LU.first, LU.second); check_results( @@ -180,9 +177,11 @@ void testExtremeValueInitialization(std::size_t number_of_blocks, double initial A[i_block][i][j] = A[0][i][j]; } - LuDecompositionPolicy lud = LuDecompositionPolicy::template Create(A); + LuDecompositionPolicy lud = + LuDecompositionPolicy::template Create(A); - auto LU = LuDecompositionPolicy::template GetLUMatrices(A, initial_value); + auto LU = LuDecompositionPolicy::template GetLUMatrices( + A, initial_value); CopyToDevice(A); CopyToDevice(LU.first); @@ -212,7 +211,8 @@ void testDiagonalMatrix(std::size_t number_of_blocks) for (std::size_t i_block = 0; i_block < number_of_blocks; ++i_block) A[i_block][i][i] = get_double(); - LuDecompositionPolicy lud = LuDecompositionPolicy::template Create(A); + LuDecompositionPolicy lud = + LuDecompositionPolicy::template Create(A); auto LU = LuDecompositionPolicy::template GetLUMatrices(A, 0); lud.template Decompose(A, LU.first, LU.second); check_results( diff --git a/test/unit/util/test_matrix.cpp b/test/unit/util/test_matrix.cpp index 096af54a9..090fff0d1 100644 --- a/test/unit/util/test_matrix.cpp +++ b/test/unit/util/test_matrix.cpp @@ -106,4 +106,14 @@ TEST(Matrix, ForEach) TEST(Matrix, SetScaler) { testSetScalar(); +} + +TEST(Matrix, Max) +{ + testMax(); +} + +TEST(Matrix, Min) +{ + testMin(); } \ No newline at end of file diff --git a/test/unit/util/test_matrix_policy.hpp b/test/unit/util/test_matrix_policy.hpp index 29bb84aea..812e61689 100644 --- a/test/unit/util/test_matrix_policy.hpp +++ b/test/unit/util/test_matrix_policy.hpp @@ -277,5 +277,57 @@ MatrixPolicy testSetScalar() EXPECT_EQ(elem, 2.0); } + return matrix; +} + +template class MatrixPolicy> +MatrixPolicy testMax() +{ + MatrixPolicy matrix{ 2, 3, 0.0 }; + + matrix.Max(2.0); + + for (auto &elem : matrix.AsVector()) + { + EXPECT_EQ(elem, 2.0); + } + + matrix = 1.0; + matrix[1][1] = 3.0; + matrix.Max(2.0); + + EXPECT_EQ(matrix[0][0], 2.0); + EXPECT_EQ(matrix[0][1], 2.0); + EXPECT_EQ(matrix[0][2], 2.0); + EXPECT_EQ(matrix[1][0], 2.0); + EXPECT_EQ(matrix[1][1], 3.0); + EXPECT_EQ(matrix[1][2], 2.0); + + return matrix; +} + +template class MatrixPolicy> +MatrixPolicy testMin() +{ + MatrixPolicy matrix{ 2, 3, 0.0 }; + + matrix.Min(2.0); + + for (auto &elem : matrix.AsVector()) + { + EXPECT_EQ(elem, 0.0); + } + + matrix = 1.0; + matrix[1][1] = 3.0; + matrix.Min(2.0); + + EXPECT_EQ(matrix[0][0], 1.0); + EXPECT_EQ(matrix[0][1], 1.0); + EXPECT_EQ(matrix[0][2], 1.0); + EXPECT_EQ(matrix[1][0], 1.0); + EXPECT_EQ(matrix[1][1], 2.0); + EXPECT_EQ(matrix[1][2], 1.0); + return matrix; } \ No newline at end of file diff --git a/test/unit/util/test_sparse_matrix_standard_ordering_column.cpp b/test/unit/util/test_sparse_matrix_standard_ordering_column.cpp index 81cd3f7e9..f5d5f11fb 100644 --- a/test/unit/util/test_sparse_matrix_standard_ordering_column.cpp +++ b/test/unit/util/test_sparse_matrix_standard_ordering_column.cpp @@ -31,7 +31,7 @@ TEST(SparseCompressedColumnMatrix, SingleBlockMatrix) { { auto matrix = testSingleBlockMatrix(); - + { std::size_t elem = matrix.VectorIndex(3, 2); EXPECT_EQ(elem, 3); diff --git a/test/unit/util/test_sparse_matrix_standard_ordering_row.cpp b/test/unit/util/test_sparse_matrix_standard_ordering_row.cpp index e809f44c3..e40307a38 100644 --- a/test/unit/util/test_sparse_matrix_standard_ordering_row.cpp +++ b/test/unit/util/test_sparse_matrix_standard_ordering_row.cpp @@ -31,7 +31,7 @@ TEST(SparseCompressedRowMatrix, SingleBlockMatrix) { { auto matrix = testSingleBlockMatrix(); - + { std::size_t elem = matrix.VectorIndex(3, 2); EXPECT_EQ(elem, 4); diff --git a/test/unit/util/test_vector_matrix.cpp b/test/unit/util/test_vector_matrix.cpp index 799eff40c..061689f09 100644 --- a/test/unit/util/test_vector_matrix.cpp +++ b/test/unit/util/test_vector_matrix.cpp @@ -105,4 +105,20 @@ TEST(VectorMatrix, SetScaler) testSetScalar(); testSetScalar(); testSetScalar(); +} + +TEST(VectorMatrix, Max) +{ + testMax(); + testMax(); + testMax(); + testMax(); +} + +TEST(VectorMatrix, Min) +{ + testMin(); + testMin(); + testMin(); + testMin(); } \ No newline at end of file