diff --git a/include/micm/cuda/solver/cuda_lu_decomposition.cuh b/include/micm/cuda/solver/cuda_lu_decomposition.cuh index b3e619b5b..22cf7cb3d 100644 --- a/include/micm/cuda/solver/cuda_lu_decomposition.cuh +++ b/include/micm/cuda/solver/cuda_lu_decomposition.cuh @@ -14,8 +14,7 @@ namespace micm const CudaMatrixParam& A_param, CudaMatrixParam& L_param, CudaMatrixParam& U_param, - const LuDecomposeParam& devstruct, - bool& is_singular); + const LuDecomposeParam& devstruct); /// This is the function that will copy the constant data /// members of class "CudaLuDecomposition" to the device; diff --git a/include/micm/cuda/solver/cuda_lu_decomposition.hpp b/include/micm/cuda/solver/cuda_lu_decomposition.hpp index c8058a14e..7acd23444 100644 --- a/include/micm/cuda/solver/cuda_lu_decomposition.hpp +++ b/include/micm/cuda/solver/cuda_lu_decomposition.hpp @@ -85,14 +85,6 @@ namespace micm /// @param A is the sparse matrix to decompose /// @param L is the lower triangular matrix created by decomposition /// @param U is the upper triangular matrix created by decomposition - /// @param is_singular Flag that is set to true if A is singular; false otherwise - template - requires(CudaMatrix&& VectorizableSparse) void Decompose( - const SparseMatrixPolicy& A, - SparseMatrixPolicy& L, - SparseMatrixPolicy& U, - bool& is_singular) const; - template requires(CudaMatrix&& VectorizableSparse) void Decompose( const SparseMatrixPolicy& A, @@ -100,27 +92,14 @@ namespace micm SparseMatrixPolicy& U) const; }; - template - requires(CudaMatrix&& VectorizableSparse) void CudaLuDecomposition::Decompose( - const SparseMatrixPolicy& A, - SparseMatrixPolicy& L, - SparseMatrixPolicy& U, - bool& is_singular) 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 - micm::cuda::DecomposeKernelDriver(A.AsDeviceParam(), L_param, U_param, this->devstruct_, is_singular); - } - template requires(CudaMatrix&& VectorizableSparse) void CudaLuDecomposition::Decompose( const SparseMatrixPolicy& A, SparseMatrixPolicy& L, SparseMatrixPolicy& U) const { - bool is_singular = false; 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 - micm::cuda::DecomposeKernelDriver(A.AsDeviceParam(), L_param, U_param, this->devstruct_, is_singular); + micm::cuda::DecomposeKernelDriver(A.AsDeviceParam(), L_param, U_param, this->devstruct_); } } // end of namespace micm diff --git a/include/micm/cuda/util/cuda_param.hpp b/include/micm/cuda/util/cuda_param.hpp index d415a3fa2..766587c14 100644 --- a/include/micm/cuda/util/cuda_param.hpp +++ b/include/micm/cuda/util/cuda_param.hpp @@ -33,7 +33,6 @@ struct ProcessSetParam /// This struct could be allocated on the host or device; struct LuDecomposeParam { - bool* is_singular = nullptr; std::pair* niLU_ = nullptr; char* do_aik_ = nullptr; std::size_t* aik_ = nullptr; diff --git a/include/micm/jit/solver/jit_linear_solver.hpp b/include/micm/jit/solver/jit_linear_solver.hpp index 2572f17c7..3691cd922 100644 --- a/include/micm/jit/solver/jit_linear_solver.hpp +++ b/include/micm/jit/solver/jit_linear_solver.hpp @@ -42,12 +42,10 @@ 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 - /// @param is_singular Flag that will be set to true if matrix is singular; false otherwise void Factor( SparseMatrixPolicy& matrix, SparseMatrixPolicy& lower_matrix, - SparseMatrixPolicy& upper_matrix, - bool& is_singular) const; + SparseMatrixPolicy& upper_matrix) const; /// @brief Solve for x in Ax = b template diff --git a/include/micm/jit/solver/jit_linear_solver.inl b/include/micm/jit/solver/jit_linear_solver.inl index 545fa7318..b678aeac1 100644 --- a/include/micm/jit/solver/jit_linear_solver.inl +++ b/include/micm/jit/solver/jit_linear_solver.inl @@ -60,10 +60,9 @@ namespace micm inline void JitLinearSolver::Factor( SparseMatrixPolicy &matrix, SparseMatrixPolicy &lower_matrix, - SparseMatrixPolicy &upper_matrix, - bool &is_singular) const + SparseMatrixPolicy &upper_matrix) const { - LinearSolver::Factor(matrix, lower_matrix, upper_matrix, is_singular); + LinearSolver::Factor(matrix, lower_matrix, upper_matrix); } template diff --git a/include/micm/jit/solver/jit_lu_decomposition.hpp b/include/micm/jit/solver/jit_lu_decomposition.hpp index 5be4795db..20f721a6a 100644 --- a/include/micm/jit/solver/jit_lu_decomposition.hpp +++ b/include/micm/jit/solver/jit_lu_decomposition.hpp @@ -39,9 +39,8 @@ namespace micm /// @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 - /// @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) + void Decompose(const SparseMatrixPolicy &A, SparseMatrixPolicy &lower, SparseMatrixPolicy &upper) const; private: diff --git a/include/micm/jit/solver/jit_lu_decomposition.inl b/include/micm/jit/solver/jit_lu_decomposition.inl index 452e8c03f..d7e3e9ca3 100644 --- a/include/micm/jit/solver/jit_lu_decomposition.inl +++ b/include/micm/jit/solver/jit_lu_decomposition.inl @@ -209,22 +209,12 @@ namespace micm void JitLuDecomposition::Decompose( const SparseMatrixPolicy &A, SparseMatrixPolicy &lower, - SparseMatrixPolicy &upper, - bool &is_singular) const + SparseMatrixPolicy &upper) const { - is_singular = false; decompose_function_(A.AsVector().data(), lower.AsVector().data(), upper.AsVector().data()); for (size_t block = 0; block < A.NumberOfBlocks(); ++block) { auto diagonals = upper.DiagonalIndices(block); - for (const auto &diagonal : diagonals) - { - if (upper.AsVector()[diagonal] == 0) - { - is_singular = true; - return; - } - } } } } // namespace micm \ No newline at end of file diff --git a/include/micm/solver/backward_euler.inl b/include/micm/solver/backward_euler.inl index ec7e4b9aa..25fe68353 100644 --- a/include/micm/solver/backward_euler.inl +++ b/include/micm/solver/backward_euler.inl @@ -66,8 +66,6 @@ namespace micm std::size_t n_successful_integrations = 0; std::size_t n_convergence_failures = 0; - bool singular = false; - auto derived_class_temporary_variables = static_cast*>(state.temporary_variables_.get()); auto& Yn = derived_class_temporary_variables->Yn_; @@ -112,7 +110,7 @@ namespace micm // (y_{n+1} - y_n) / H = f(t_{n+1}, y_{n+1}) // try to find the root by factoring and solving the linear system - linear_solver_.Factor(state.jacobian_, state.lower_matrix_, state.upper_matrix_, singular); + linear_solver_.Factor(state.jacobian_, state.lower_matrix_, state.upper_matrix_); result.stats_.decompositions_++; // forcing_blk in camchem diff --git a/include/micm/solver/linear_solver.hpp b/include/micm/solver/linear_solver.hpp index d028546c6..1f72962c3 100644 --- a/include/micm/solver/linear_solver.hpp +++ b/include/micm/solver/linear_solver.hpp @@ -78,12 +78,10 @@ namespace micm /// @brief Decompose the matrix into upper and lower triangular matrices /// @param matrix Matrix to decompose into lower and upper triangular matrices - /// @param is_singular Flag that is set to true if matrix is singular; false otherwise void Factor( const SparseMatrixPolicy& matrix, SparseMatrixPolicy& lower_matrix, - SparseMatrixPolicy& upper_matrix, - bool& is_singular) const; + SparseMatrixPolicy& upper_matrix) const; /// @brief Solve for x in Ax = b. x should be a copy of b and after Solve finishes x will contain the result template diff --git a/include/micm/solver/linear_solver.inl b/include/micm/solver/linear_solver.inl index f1b2a8492..a04947287 100644 --- a/include/micm/solver/linear_solver.inl +++ b/include/micm/solver/linear_solver.inl @@ -111,12 +111,11 @@ namespace micm inline void LinearSolver::Factor( const SparseMatrixPolicy& matrix, SparseMatrixPolicy& lower_matrix, - SparseMatrixPolicy& upper_matrix, - bool& is_singular) const + SparseMatrixPolicy& upper_matrix) const { MICM_PROFILE_FUNCTION(); - lu_decomp_.template Decompose(matrix, lower_matrix, upper_matrix, is_singular); + lu_decomp_.template Decompose(matrix, lower_matrix, upper_matrix); } template diff --git a/include/micm/solver/lu_decomposition.hpp b/include/micm/solver/lu_decomposition.hpp index 0b478a4ec..42b8061c9 100644 --- a/include/micm/solver/lu_decomposition.hpp +++ b/include/micm/solver/lu_decomposition.hpp @@ -122,19 +122,16 @@ namespace micm /// @param A Sparse matrix to decompose /// @param L The lower triangular matrix created by decomposition /// @param U The upper triangular matrix created by decomposition - /// @param is_singular Flag that is set to true if A is singular; false otherwise template requires(!VectorizableSparse) void Decompose( const SparseMatrixPolicy& A, SparseMatrixPolicy& L, - SparseMatrixPolicy& U, - bool& is_singular) const; + SparseMatrixPolicy& U) const; template requires(VectorizableSparse) void Decompose( const SparseMatrixPolicy& A, SparseMatrixPolicy& L, - SparseMatrixPolicy& U, - bool& is_singular) const; + SparseMatrixPolicy& U) const; private: /// @brief Initialize arrays for the LU decomposition diff --git a/include/micm/solver/lu_decomposition.inl b/include/micm/solver/lu_decomposition.inl index 6385e53e0..54137fe5f 100644 --- a/include/micm/solver/lu_decomposition.inl +++ b/include/micm/solver/lu_decomposition.inl @@ -168,11 +168,9 @@ namespace micm requires(!VectorizableSparse) inline void LuDecomposition::Decompose( const SparseMatrixPolicy& A, SparseMatrixPolicy& L, - SparseMatrixPolicy& U, - bool& is_singular) const + SparseMatrixPolicy& U) const { MICM_PROFILE_FUNCTION(); - is_singular = false; // Loop over blocks for (std::size_t i_block = 0; i_block < A.NumberOfBlocks(); ++i_block) @@ -226,21 +224,11 @@ namespace micm L_vector[lki_nkj->first] -= L_vector[lkj_uji->first] * U_vector[lkj_uji->second]; ++lkj_uji; } - - if (U_vector[*uii] == 0.0) - { - is_singular = true; - } L_vector[lki_nkj->first] /= U_vector[*uii]; ++lki_nkj; ++uii; } } - // check the bottom right corner of the matrix - if (U_vector[*uii] == 0.0) - { - is_singular = true; - } } } @@ -248,8 +236,7 @@ namespace micm requires(VectorizableSparse) inline void LuDecomposition::Decompose( const SparseMatrixPolicy& A, SparseMatrixPolicy& L, - SparseMatrixPolicy& U, - bool& is_singular) const + SparseMatrixPolicy& U) const { MICM_PROFILE_FUNCTION(); @@ -258,7 +245,6 @@ namespace micm const std::size_t A_GroupSizeOfFlatBlockSize = A.GroupSize(A.FlatBlockSize()); const std::size_t L_GroupSizeOfFlatBlockSize = L.GroupSize(L.FlatBlockSize()); const std::size_t U_GroupSizeOfFlatBlockSize = U.GroupSize(U.FlatBlockSize()); - is_singular = false; // Loop over groups of blocks for (std::size_t i_group = 0; i_group < A.NumberOfGroups(A_BlockSize); ++i_group) @@ -328,29 +314,12 @@ namespace micm const std::size_t uii_deref = *uii; for (std::size_t i_cell = 0; i_cell < n_cells; ++i_cell) { - if (U_vector[uii_deref + i_cell] == 0.0) - { - is_singular = true; - } L_vector[lki_nkj_first + i_cell] /= U_vector[uii_deref + i_cell]; } ++lki_nkj; ++uii; } } - const std::size_t uii_deref = *uii; - if (n_cells != A_GroupVectorSize) - { - // check the bottom right corner of the matrix - for (std::size_t i_cell = 0; i_cell < n_cells; ++i_cell) - { - if (U_vector[uii_deref + i_cell] == 0.0) - { - is_singular = true; - break; - } - } - } } } diff --git a/include/micm/solver/rosenbrock.hpp b/include/micm/solver/rosenbrock.hpp index 714dce978..0955a3b32 100644 --- a/include/micm/solver/rosenbrock.hpp +++ b/include/micm/solver/rosenbrock.hpp @@ -105,14 +105,12 @@ namespace micm /// @brief Perform the LU decomposition of the matrix /// @param H The time step /// @param gamma The gamma value - /// @param singular A flag to indicate if the matrix is singular /// @param number_densities The number densities /// @param stats The solver stats /// @param state The state void LinearFactor( double& H, const double gamma, - bool& singular, const auto& number_densities, SolverStats& stats, auto& state) const; diff --git a/include/micm/solver/rosenbrock.inl b/include/micm/solver/rosenbrock.inl index 3b4127a57..bb2226819 100644 --- a/include/micm/solver/rosenbrock.inl +++ b/include/micm/solver/rosenbrock.inl @@ -69,14 +69,8 @@ namespace micm // Repeat step calculation until current step accepted while (!accepted) { - bool is_singular{ false }; // Form and factor the rosenbrock ode jacobian - LinearFactor(H, parameters_.gamma_[0], is_singular, Y, stats, state); - if (is_singular) - { - result.state_ = SolverState::RepeatedlySingularMatrix; - break; - } + LinearFactor(H, parameters_.gamma_[0], Y, stats, state); // Compute the stages for (uint64_t stage = 0; stage < parameters_.stages_; ++stage) @@ -230,37 +224,17 @@ namespace micm inline void AbstractRosenbrockSolver::LinearFactor( double& H, const double gamma, - bool& singular, const auto& number_densities, SolverStats& stats, auto& state) const { MICM_PROFILE_FUNCTION(); - uint64_t n_consecutive = 0; - singular = false; - while (true) - { - double alpha = 1 / (H * gamma); - static_cast(this)->AlphaMinusJacobian(state.jacobian_, alpha); - - linear_solver_.Factor(state.jacobian_, state.lower_matrix_, state.upper_matrix_, singular); - stats.decompositions_ += 1; + double alpha = 1 / (H * gamma); + static_cast(this)->AlphaMinusJacobian(state.jacobian_, alpha); - // if we are checking for singularity and the matrix is not singular, we can break the loop - // if we are not checking for singularity, we always break the loop - if (!singular || !parameters_.check_singularity_) - break; - - stats.singular_ += 1; - if (++n_consecutive > 5) - break; - H /= 2; - // Reconstruct the Jacobian matrix if a substepping is performed here - state.jacobian_.Fill(0); - rates_.SubtractJacobianTerms(state.rate_constants_, number_densities, state.jacobian_); - stats.jacobian_updates_ += 1; - } + linear_solver_.Factor(state.jacobian_, state.lower_matrix_, state.upper_matrix_); + stats.decompositions_ += 1; } template diff --git a/include/micm/solver/rosenbrock_solver_parameters.hpp b/include/micm/solver/rosenbrock_solver_parameters.hpp index 6ae0a6b6f..1c5e6037c 100644 --- a/include/micm/solver/rosenbrock_solver_parameters.hpp +++ b/include/micm/solver/rosenbrock_solver_parameters.hpp @@ -63,8 +63,6 @@ namespace micm std::vector absolute_tolerance_; double relative_tolerance_{ 1e-6 }; - bool check_singularity_{ false }; // Check for singular A matrix in linear solve of A x = b - // Print RosenbrockSolverParameters to console void Print() const; diff --git a/include/micm/solver/solver_result.hpp b/include/micm/solver/solver_result.hpp index e9bedb23b..fb9797df5 100644 --- a/include/micm/solver/solver_result.hpp +++ b/include/micm/solver/solver_result.hpp @@ -45,9 +45,6 @@ namespace micm uint64_t decompositions_{}; /// @brief The number of linear solves uint64_t solves_{}; - /// @brief The number of times a singular matrix is detected. For now, this will always be zero as we assume the matrix - /// is never singular - uint64_t singular_{}; /// @brief Set all member variables to zero void Reset(); @@ -62,7 +59,6 @@ namespace micm rejected_ = 0; decompositions_ = 0; solves_ = 0; - singular_ = 0; } inline std::string SolverStateToString(const SolverState& state) diff --git a/src/process/process_set.cu b/src/process/process_set.cu index 67c9e3431..7c165934d 100644 --- a/src/process/process_set.cu +++ b/src/process/process_set.cu @@ -19,17 +19,17 @@ namespace micm // Local device variables size_t react_id_offset, prod_id_offset, yield_offset; - size_t* d_number_of_reactants = devstruct.number_of_reactants_; - size_t* d_reactant_ids = devstruct.reactant_ids_; - size_t* d_number_of_products = devstruct.number_of_products_; - size_t* d_product_ids = devstruct.product_ids_; - double* d_yields = devstruct.yields_; + const size_t* const d_number_of_reactants = devstruct.number_of_reactants_; + const size_t* const d_reactant_ids = devstruct.reactant_ids_; + const size_t* const d_number_of_products = devstruct.number_of_products_; + const size_t* const d_product_ids = devstruct.product_ids_; + const double* const d_yields = devstruct.yields_; const size_t number_of_grid_cells = rate_constants_param.number_of_grid_cells_; const size_t number_of_reactions = - rate_constants_param.number_of_elements_ / rate_constants_param.number_of_grid_cells_; - const double* d_rate_constants = rate_constants_param.d_data_; - const double* d_state_variables = state_variables_param.d_data_; - double* d_forcing = forcing_param.d_data_; + rate_constants_param.number_of_elements_ / rate_constants_param.number_of_grid_cells_; + const double* const d_rate_constants = rate_constants_param.d_data_; + const double* const d_state_variables = state_variables_param.d_data_; + double* const d_forcing = forcing_param.d_data_; if (tid < number_of_grid_cells) { @@ -73,17 +73,17 @@ namespace micm size_t react_ids_offset = 0; size_t yields_offset = 0; size_t flat_id_offset = 0; - size_t* d_number_of_reactants = devstruct.number_of_reactants_; - size_t* d_reactant_ids = devstruct.reactant_ids_; - size_t* d_number_of_products = devstruct.number_of_products_; - size_t* d_jacobian_flat_ids = devstruct.jacobian_flat_ids_; - double* d_yields = devstruct.yields_; + const size_t* const d_number_of_reactants = devstruct.number_of_reactants_; + const size_t* const d_reactant_ids = devstruct.reactant_ids_; + const size_t* const d_number_of_products = devstruct.number_of_products_; + const size_t* const d_jacobian_flat_ids = devstruct.jacobian_flat_ids_; + const double* const d_yields = devstruct.yields_; const size_t number_of_grid_cells = rate_constants_param.number_of_grid_cells_; const size_t number_of_reactions = rate_constants_param.number_of_elements_ / rate_constants_param.number_of_grid_cells_; - const double* d_rate_constants = rate_constants_param.d_data_; - const double* d_state_variables = state_variables_param.d_data_; - double* d_jacobian = jacobian_param.d_data_; + const double* const d_rate_constants = rate_constants_param.d_data_; + const double* const d_state_variables = state_variables_param.d_data_; + double* const d_jacobian = jacobian_param.d_data_; if (tid < number_of_grid_cells) { diff --git a/src/solver/linear_solver.cu b/src/solver/linear_solver.cu index 1ffa1ef1c..45e2f5a78 100644 --- a/src/solver/linear_solver.cu +++ b/src/solver/linear_solver.cu @@ -20,17 +20,17 @@ namespace micm size_t tid = blockIdx.x * BLOCK_SIZE + threadIdx.x; // Local device variables - std::pair* d_nLij_Lii = devstruct.nLij_Lii_; - std::pair* d_Lij_yj = devstruct.Lij_yj_; - std::pair* d_nUij_Uii = devstruct.nUij_Uii_; - std::pair* d_Uij_xj = devstruct.Uij_xj_; - size_t nLij_Lii_size = devstruct.nLij_Lii_size_; - size_t nUij_Uii_size = devstruct.nUij_Uii_size_; - - double* d_L = L_param.d_data_; - double* d_U = U_param.d_data_; - double* d_x = x_param.d_data_; - double* d_y = d_x; // Alias d_x for consistency with equation, but to reuse memory + const std::pair* const d_nLij_Lii = devstruct.nLij_Lii_; + const std::pair* const d_Lij_yj = devstruct.Lij_yj_; + const std::pair* const d_nUij_Uii = devstruct.nUij_Uii_; + const std::pair* const d_Uij_xj = devstruct.Uij_xj_; + const size_t nLij_Lii_size = devstruct.nLij_Lii_size_; + const size_t nUij_Uii_size = devstruct.nUij_Uii_size_; + + const double* const d_L = L_param.d_data_; + const double* const d_U = U_param.d_data_; + double* const d_x = x_param.d_data_; + double* const d_y = d_x; // Alias d_x for consistency with equation, but to reuse memory const size_t number_of_grid_cells = x_param.number_of_grid_cells_; const size_t number_of_species = x_param.number_of_elements_ / number_of_grid_cells; diff --git a/src/solver/lu_decomposition.cu b/src/solver/lu_decomposition.cu index aef3b957d..414b7a87f 100644 --- a/src/solver/lu_decomposition.cu +++ b/src/solver/lu_decomposition.cu @@ -18,17 +18,17 @@ namespace micm size_t tid = blockIdx.x * BLOCK_SIZE + threadIdx.x; // Local device variables - std::pair* d_niLU = devstruct.niLU_; - char* d_do_aik = devstruct.do_aik_; - size_t* d_aik = devstruct.aik_; - std::pair* d_uik_nkj = devstruct.uik_nkj_; - std::pair* d_lij_ujk = devstruct.lij_ujk_; - char* d_do_aki = devstruct.do_aki_; - size_t* d_aki = devstruct.aki_; - std::pair* d_lki_nkj = devstruct.lki_nkj_; - std::pair* d_lkj_uji = devstruct.lkj_uji_; - size_t* d_uii = devstruct.uii_; - size_t niLU_size = devstruct.niLU_size_; + const std::pair* const d_niLU = devstruct.niLU_; + const char* const d_do_aik = devstruct.do_aik_; + const size_t* const d_aik = devstruct.aik_; + const std::pair* const d_uik_nkj = devstruct.uik_nkj_; + const std::pair* const d_lij_ujk = devstruct.lij_ujk_; + const char* const d_do_aki = devstruct.do_aki_; + const size_t* const d_aki = devstruct.aki_; + const std::pair* const d_lki_nkj = devstruct.lki_nkj_; + const std::pair* const d_lkj_uji = devstruct.lkj_uji_; + const size_t* const d_uii = devstruct.uii_; + const size_t niLU_size = devstruct.niLU_size_; size_t do_aik_offset = 0; size_t aik_offset = 0; @@ -40,12 +40,10 @@ namespace micm size_t lki_nkj_offset = 0; size_t uii_offset = 0; - double* d_A = A_param.d_data_; - double* d_L = L_param.d_data_; - double* d_U = U_param.d_data_; - size_t number_of_grid_cells = A_param.number_of_grid_cells_; - bool* d_is_singular = devstruct.is_singular; - *d_is_singular = false; + const double* const d_A = A_param.d_data_; + double* const d_L = L_param.d_data_; + double* const d_U = U_param.d_data_; + const size_t number_of_grid_cells = A_param.number_of_grid_cells_; if (tid < number_of_grid_cells) { @@ -101,20 +99,11 @@ namespace micm d_L[L_idx_1] -= d_L[L_idx_2] * d_U[U_idx]; ++lkj_uji_offset; } - if (d_U[d_uii[uii_offset] + tid] == 0.0) - { - *d_is_singular = true; - } d_L[d_lki_nkj[lki_nkj_offset].first + tid] /= d_U[d_uii[uii_offset] + tid]; ++lki_nkj_offset; ++uii_offset; } } - // check the bottom right corner of the matrix - if (d_U[d_uii[uii_offset] + tid] == 0.0) - { - *d_is_singular = true; - } } } // end of CUDA kernel @@ -172,10 +161,6 @@ namespace micm CHECK_CUDA_ERROR( cudaMallocAsync(&(devstruct.uii_), uii_bytes, micm::cuda::CudaStreamSingleton::GetInstance().GetCudaStream(0)), "cudaMalloc"); - CHECK_CUDA_ERROR( - cudaMallocAsync( - &devstruct.is_singular, sizeof(bool), micm::cuda::CudaStreamSingleton::GetInstance().GetCudaStream(0)), - "cudaMalloc"); /// Copy the data from host to device CHECK_CUDA_ERROR( @@ -267,10 +252,6 @@ namespace micm /// members of class "CudaLuDecomposition" on the device void FreeConstData(LuDecomposeParam& devstruct) { - if (devstruct.is_singular != nullptr) - CHECK_CUDA_ERROR( - cudaFreeAsync(devstruct.is_singular, micm::cuda::CudaStreamSingleton::GetInstance().GetCudaStream(0)), - "cudaFree"); if (devstruct.niLU_ != nullptr) CHECK_CUDA_ERROR( cudaFreeAsync(devstruct.niLU_, micm::cuda::CudaStreamSingleton::GetInstance().GetCudaStream(0)), "cudaFree"); @@ -307,20 +288,12 @@ namespace micm const CudaMatrixParam& A_param, CudaMatrixParam& L_param, CudaMatrixParam& U_param, - const LuDecomposeParam& devstruct, - bool& is_singular) + const LuDecomposeParam& devstruct) { // Launch the CUDA kernel for LU decomposition size_t number_of_blocks = (A_param.number_of_grid_cells_ + BLOCK_SIZE - 1) / BLOCK_SIZE; DecomposeKernel<<>>( A_param, L_param, U_param, devstruct); - // Copy the boolean result from device back to host - cudaMemcpyAsync( - &is_singular, - devstruct.is_singular, - sizeof(bool), - cudaMemcpyDeviceToHost, - micm::cuda::CudaStreamSingleton::GetInstance().GetCudaStream(0)); } // end of DecomposeKernelDriver } // end of namespace cuda } // end of namespace micm diff --git a/src/solver/rosenbrock.cu b/src/solver/rosenbrock.cu index 59d59f46b..6a6ab6961 100644 --- a/src/solver/rosenbrock.cu +++ b/src/solver/rosenbrock.cu @@ -142,11 +142,11 @@ namespace micm const size_t n, bool is_first_call) { - double* d_y_old = y_old_param.d_data_; - double* d_y_new = y_new_param.d_data_; - double* d_errors_input = devstruct.errors_input_; - double* d_errors_output = devstruct.errors_output_; - const double* atol = devstruct.absolute_tolerance_; + const double* const d_y_old = y_old_param.d_data_; + const double* const d_y_new = y_new_param.d_data_; + double* const d_errors_input = devstruct.errors_input_; + double* const d_errors_output = devstruct.errors_output_; + const double* const atol = devstruct.absolute_tolerance_; const double rtol = ros_param.relative_tolerance_; const size_t number_of_grid_cells = y_old_param.number_of_grid_cells_; @@ -250,11 +250,11 @@ namespace micm { // Local device variables double d_ymax, d_scale; - double* d_y_old = y_old_param.d_data_; - double* d_y_new = y_new_param.d_data_; - double* d_errors = devstruct.errors_input_; - const double* atol = devstruct.absolute_tolerance_; - double rtol = ros_param.relative_tolerance_; + const double* const d_y_old = y_old_param.d_data_; + const double* const d_y_new = y_new_param.d_data_; + double* const d_errors = devstruct.errors_input_; + const double* const atol = devstruct.absolute_tolerance_; + const double rtol = ros_param.relative_tolerance_; const size_t num_elements = devstruct.errors_size_; const size_t number_of_grid_cells = y_old_param.number_of_grid_cells_; diff --git a/test/integration/test_terminator.cpp b/test/integration/test_terminator.cpp index b86d35cb5..0384301de 100644 --- a/test/integration/test_terminator.cpp +++ b/test/integration/test_terminator.cpp @@ -32,7 +32,6 @@ TEST(RosenbrockSolver, Terminator) auto builder = micm::CpuSolverBuilder(parameters).SetIgnoreUnusedSpecies(true); TestTerminator(builder, 4); } - parameters.check_singularity_ = true; { auto builder = micm::CpuSolverBuilder(parameters).SetIgnoreUnusedSpecies(true); TestTerminator(builder, 1); diff --git a/test/regression/RosenbrockChapman/chapman_ode_solver.hpp b/test/regression/RosenbrockChapman/chapman_ode_solver.hpp index c565829d0..16d4c5d16 100644 --- a/test/regression/RosenbrockChapman/chapman_ode_solver.hpp +++ b/test/regression/RosenbrockChapman/chapman_ode_solver.hpp @@ -72,7 +72,6 @@ namespace micm uint64_t rejected{}; // Nrej uint64_t decompositions{}; // Ndec uint64_t solves{}; // Nsol - uint64_t singular{}; // Nsng uint64_t total_steps{}; // Ntotstp void Reset(); @@ -173,11 +172,9 @@ namespace micm /// @param H time step (seconds) /// @param gamma time step factor for specific rosenbrock method /// @param Y constituent concentration (molec/cm^3) - /// @param singular indicates if the matrix is singular std::vector lin_factor( double& H, const double& gamma, - bool& singular, const std::vector& number_densities, const double& number_density_air, const std::vector& rate_constants); @@ -209,7 +206,6 @@ namespace micm rejected = 0; decompositions = 0; solves = 0; - singular = 0; total_steps = 0; } @@ -365,15 +361,9 @@ namespace micm { break; } - bool is_singular{ false }; // Form and factor the rosenbrock ode jacobian - auto ode_jacobian = lin_factor(H, parameters_.gamma_[0], is_singular, Y, number_density_air, rate_constants); + auto ode_jacobian = lin_factor(H, parameters_.gamma_[0], Y, number_density_air, rate_constants); stats_.jacobian_updates += 1; - if (is_singular) - { - result.state_ = SolverState::RepeatedlySingularMatrix; - break; - } // Compute the stages for (uint64_t stage = 0; stage < parameters_.stages_; ++stage) @@ -787,7 +777,6 @@ namespace micm inline std::vector ChapmanODESolver::lin_factor( double& H, const double& gamma, - bool& singular, const std::vector& number_densities, const double& number_density_air, const std::vector& rate_constants) @@ -801,35 +790,11 @@ namespace micm std::function)> is_successful = [](const std::vector& jacobian) { return true; }; std::vector ode_jacobian; uint64_t n_consecutive = 0; - singular = true; - - while (true) - { - double alpha = 1 / (H * gamma); - // compute jacobian decomposition of alpha*I - dforce_dy - ode_jacobian = factored_alpha_minus_jac(dforce_dy(rate_constants, number_densities, number_density_air), alpha); - stats_.decompositions += 1; - - if (is_successful(ode_jacobian)) - { - singular = false; - break; - } - else - { - stats_.singular += 1; - n_consecutive += 1; - if (n_consecutive <= 5) - { - H /= 2; - } - else - { - break; - } - } - } + double alpha = 1 / (H * gamma); + // compute jacobian decomposition of alpha*I - dforce_dy + ode_jacobian = factored_alpha_minus_jac(dforce_dy(rate_constants, number_densities, number_density_air), alpha); + stats_.decompositions += 1; return ode_jacobian; } diff --git a/test/tutorial/test_jit_tutorial.cpp b/test/tutorial/test_jit_tutorial.cpp index 2faad56ff..2609e0bea 100644 --- a/test/tutorial/test_jit_tutorial.cpp +++ b/test/tutorial/test_jit_tutorial.cpp @@ -61,7 +61,6 @@ auto run_solver(auto& solver) total_stats.rejected_ += result.stats_.rejected_; total_stats.decompositions_ += result.stats_.decompositions_; total_stats.solves_ += result.stats_.solves_; - total_stats.singular_ += result.stats_.singular_; } } @@ -132,7 +131,6 @@ int main(const int argc, const char* argv[]) std::cout << "\trejected: " << result_stats.rejected_ << std::endl; std::cout << "\tdecompositions: " << result_stats.decompositions_ << std::endl; std::cout << "\tsolves: " << result_stats.solves_ << std::endl; - std::cout << "\tsingular: " << result_stats.singular_ << std::endl; auto jit_result_stats = std::get<1>(jit_result_tuple); std::cout << "JIT solve stats: " << std::endl; @@ -144,7 +142,6 @@ int main(const int argc, const char* argv[]) std::cout << "\trejected: " << jit_result_stats.rejected_ << std::endl; std::cout << "\tdecompositions: " << jit_result_stats.decompositions_ << std::endl; std::cout << "\tsolves: " << jit_result_stats.solves_ << std::endl; - std::cout << "\tsingular: " << jit_result_stats.singular_ << std::endl; auto result = std::get<0>(result_tuple); auto jit_result = std::get<0>(jit_result_tuple); diff --git a/test/tutorial/test_solver_configuration.cpp b/test/tutorial/test_solver_configuration.cpp index f13125c25..522908d23 100644 --- a/test/tutorial/test_solver_configuration.cpp +++ b/test/tutorial/test_solver_configuration.cpp @@ -64,7 +64,6 @@ void test_solver_type(auto& solver) total_stats.rejected_ += result.stats_.rejected_; total_stats.decompositions_ += result.stats_.decompositions_; total_stats.solves_ += result.stats_.solves_; - total_stats.singular_ += result.stats_.singular_; elapsed_solve_time = result.final_time_; } @@ -80,7 +79,6 @@ void test_solver_type(auto& solver) std::cout << "rejected: " << total_stats.rejected_ << std::endl; std::cout << "decompositions: " << total_stats.decompositions_ << std::endl; std::cout << "solves: " << total_stats.solves_ << std::endl; - std::cout << "singular: " << total_stats.singular_ << std::endl; } int main() diff --git a/test/tutorial/test_solver_results.cpp b/test/tutorial/test_solver_results.cpp index 5c189de67..c953fc00e 100644 --- a/test/tutorial/test_solver_results.cpp +++ b/test/tutorial/test_solver_results.cpp @@ -82,6 +82,5 @@ int main() std::cout << "rejected: " << result.stats_.rejected_ << std::endl; std::cout << "decompositions: " << result.stats_.decompositions_ << std::endl; std::cout << "solves: " << result.stats_.solves_ << std::endl; - std::cout << "singular: " << result.stats_.singular_ << std::endl; std::cout << "final simulation time: " << result.final_time_ << std::endl; } \ No newline at end of file diff --git a/test/unit/cuda/solver/test_cuda_linear_solver.cpp b/test/unit/cuda/solver/test_cuda_linear_solver.cpp index d2c714832..3a7e75932 100644 --- a/test/unit/cuda/solver/test_cuda_linear_solver.cpp +++ b/test/unit/cuda/solver/test_cuda_linear_solver.cpp @@ -75,8 +75,7 @@ std::vector linearSolverGenerator(std::size_t number_of_blocks) CopyToDeviceSparse(lower_matrix); CopyToDeviceSparse(upper_matrix); - bool singular{ false }; - solver.Factor(A, lower_matrix, upper_matrix, singular); + solver.Factor(A, lower_matrix, upper_matrix); solver.template Solve(x, lower_matrix, upper_matrix); // Only copy the data to the host when it is a CudaMatrix diff --git a/test/unit/cuda/solver/test_cuda_lu_decomposition.cpp b/test/unit/cuda/solver/test_cuda_lu_decomposition.cpp index 132cd30ab..037f7c22f 100644 --- a/test/unit/cuda/solver/test_cuda_lu_decomposition.cpp +++ b/test/unit/cuda/solver/test_cuda_lu_decomposition.cpp @@ -59,8 +59,7 @@ void testCudaRandomMatrix(size_t n_grids) micm::LuDecomposition cpu_lud = micm::LuDecomposition::Create(cpu_A); auto cpu_LU = micm::LuDecomposition::GetLUMatrices(cpu_A, 0); - bool singular{ false }; - cpu_lud.Decompose(cpu_A, cpu_LU.first, cpu_LU.second, singular); + cpu_lud.Decompose(cpu_A, cpu_LU.first, cpu_LU.second); // checking GPU result again CPU size_t L_size = cpu_LU.first.AsVector().size(); diff --git a/test/unit/cuda/solver/test_cuda_rosenbrock.cpp b/test/unit/cuda/solver/test_cuda_rosenbrock.cpp index f40d11ed5..73586f7e1 100644 --- a/test/unit/cuda/solver/test_cuda_rosenbrock.cpp +++ b/test/unit/cuda/solver/test_cuda_rosenbrock.cpp @@ -251,84 +251,4 @@ TEST(RosenbrockSolver, CudaNormalizedError) testNormalizedErrorDiff<200041>(); testNormalizedErrorDiff<421875>(); testNormalizedErrorDiff<3395043>(); -} - -TEST(RosenbrockSolver, SingularSystemZeroInBottomRightOfU) -{ - auto params = micm::CudaRosenbrockSolverParameters::ThreeStageRosenbrockParameters(); - params.check_singularity_ = true; - auto vector = GpuBuilder<4>(params); - - auto vector_solver = getSingularSystemZeroInBottomRightOfU(vector).SetNumberOfGridCells(4).Build(); - - auto vector_state = vector_solver.GetState(); - - double k1 = -2; - double k2 = 1.0; - - vector_state.SetCustomRateParameter("r1", { k1, k1, k1, k1 }); - vector_state.SetCustomRateParameter("r2", { k2, k2, k2, k2 }); - - vector_state.variables_[0] = { 1.0, 1.0 }; - vector_state.variables_[1] = { 1.0, 1.0 }; - vector_state.variables_[2] = { 1.0, 1.0 }; - vector_state.variables_[3] = { 1.0, 1.0 }; - - // to get a jacobian with an LU factorization that contains a zero on the diagonal - // of U, we need det(alpha * I - jacobian) = 0 - // for the system above, that means we have to have alpha + k1 + k2 = 0 - // in this case, one of the reaction rates will be negative but it's good enough to - // test the singularity check - // alpha is 1 / (H * gamma), where H is the time step and gamma is the gamma value from - // the rosenbrock paramters - // so H needs to be 1 / ( (-k1 - k2) * gamma) - // since H is positive we need -k1 -k2 to be positive, hence the smaller, negative value for k1 - double H = 1 / ((-k1 - k2) * params.gamma_[0]); - vector_solver.solver_.parameters_.h_start_ = H; - - vector_solver.CalculateRateConstants(vector_state); - vector_state.SyncInputsToDevice(); - - auto vector_result = vector_solver.Solve(2 * H, vector_state); - vector_state.SyncOutputsToHost(); - EXPECT_NE(vector_result.stats_.singular_, 0); -} - -TEST(RosenbrockSolver, SingularSystemZeroAlongDiagonalNotBottomRight) -{ - auto params = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters(); - - double k1 = -1.0; - double k2 = -1.0; - double k3 = 1.0; - - // to get a jacobian with an LU factorization that contains a zero on the diagonal - // of U, we need det(alpha * I - jacobian) = 0 - // for the system above, that means we have to set alpha = -k1, or alpha=-k2, or alpha=k3 - double H = 1 / (-k1 * params.gamma_[0]); - - params.check_singularity_ = true; - params.h_start_ = H; - - auto vector = GpuBuilder<4>(params); - - auto vector_solver = getSolverForSingularSystemOnDiagonal(vector).SetNumberOfGridCells(4).Build(); - - auto vector_state = vector_solver.GetState(); - - vector_state.SetCustomRateParameter("r1", { k1, k1, k1, k1 }); - vector_state.SetCustomRateParameter("r2", { k2, k2, k2, k2 }); - vector_state.SetCustomRateParameter("r3", { k3, k3, k3, k3 }); - - vector_state.variables_[0] = { 1.0, 1.0, 1.0 }; - vector_state.variables_[1] = { 1.0, 1.0, 1.0 }; - vector_state.variables_[2] = { 1.0, 1.0, 1.0 }; - vector_state.variables_[3] = { 1.0, 1.0, 1.0 }; - - vector_solver.CalculateRateConstants(vector_state); - vector_state.SyncInputsToDevice(); - - auto vector_result = vector_solver.Solve(2 * H, vector_state); - vector_state.SyncOutputsToHost(); - EXPECT_NE(vector_result.stats_.singular_, 0); } \ No newline at end of file diff --git a/test/unit/jit/solver/test_jit_lu_decomposition.cpp b/test/unit/jit/solver/test_jit_lu_decomposition.cpp index 531c04ab2..fbcfea9d6 100644 --- a/test/unit/jit/solver/test_jit_lu_decomposition.cpp +++ b/test/unit/jit/solver/test_jit_lu_decomposition.cpp @@ -19,14 +19,6 @@ TEST(JitLuDecomposition, DenseMatrixVectorOrdering) testDenseMatrix>(); } -TEST(JitLuDecomposition, SingularMatrixVectorOrdering) -{ - testSingularMatrix>(); - testSingularMatrix>(); - testSingularMatrix>(); - testSingularMatrix>(); -} - TEST(JitLuDecomposition, RandomMatrixVectorOrdering) { testRandomMatrix>(1); diff --git a/test/unit/jit/solver/test_jit_rosenbrock.cpp b/test/unit/jit/solver/test_jit_rosenbrock.cpp index 542bde6f3..c980246f9 100644 --- a/test/unit/jit/solver/test_jit_rosenbrock.cpp +++ b/test/unit/jit/solver/test_jit_rosenbrock.cpp @@ -72,80 +72,4 @@ TEST(JitRosenbrockSolver, MultipleInstances) run_solver(solver1); run_solver(solver2); run_solver(solver3); -} - -TEST(JitRosenbrockSolver, SingularSystemZeroInBottomRightOfU) -{ - auto params = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters(); - params.check_singularity_ = true; - auto vector = JitBuilder<4>(params); - - auto vector_solver = getSingularSystemZeroInBottomRightOfU(vector).SetNumberOfGridCells(4).Build(); - - auto vector_state = vector_solver.GetState(); - - double k1 = -2; - double k2 = 1.0; - - vector_state.SetCustomRateParameter("r1", { k1, k1, k1, k1 }); - vector_state.SetCustomRateParameter("r2", { k2, k2, k2, k2 }); - - vector_state.variables_[0] = { 1.0, 1.0 }; - vector_state.variables_[1] = { 1.0, 1.0 }; - vector_state.variables_[2] = { 1.0, 1.0 }; - vector_state.variables_[3] = { 1.0, 1.0 }; - - // to get a jacobian with an LU factorization that contains a zero on the diagonal - // of U, we need det(alpha * I - jacobian) = 0 - // for the system above, that means we have to have alpha + k1 + k2 = 0 - // in this case, one of the reaction rates will be negative but it's good enough to - // test the singularity check - // alpha is 1 / (H * gamma), where H is the time step and gamma is the gamma value from - // the rosenbrock paramters - // so H needs to be 1 / ( (-k1 - k2) * gamma) - // since H is positive we need -k1 -k2 to be positive, hence the smaller, negative value for k1 - double H = 1 / ((-k1 - k2) * params.gamma_[0]); - vector_solver.solver_.parameters_.h_start_ = H; - - vector_solver.CalculateRateConstants(vector_state); - - auto vector_result = vector_solver.Solve(2 * H, vector_state); - EXPECT_NE(vector_result.stats_.singular_, 0); -} - -TEST(JitRosenbrockSolver, SingularSystemZeroAlongDiagonalNotBottomRight) -{ - auto params = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters(); - - double k1 = -1.0; - double k2 = -1.0; - double k3 = 1.0; - - // to get a jacobian with an LU factorization that contains a zero on the diagonal - // of U, we need det(alpha * I - jacobian) = 0 - // for the system above, that means we have to set alpha = -k1, or alpha=-k2, or alpha=k3 - double H = 1 / (-k1 * params.gamma_[0]); - - params.check_singularity_ = true; - params.h_start_ = H; - - auto vector = JitBuilder<4>(params); - - auto vector_solver = getSolverForSingularSystemOnDiagonal(vector).SetNumberOfGridCells(4).Build(); - - auto vector_state = vector_solver.GetState(); - - vector_state.SetCustomRateParameter("r1", { k1, k1, k1, k1 }); - vector_state.SetCustomRateParameter("r2", { k2, k2, k2, k2 }); - vector_state.SetCustomRateParameter("r3", { k3, k3, k3, k3 }); - - vector_state.variables_[0] = { 1.0, 1.0, 1.0 }; - vector_state.variables_[1] = { 1.0, 1.0, 1.0 }; - vector_state.variables_[2] = { 1.0, 1.0, 1.0 }; - vector_state.variables_[3] = { 1.0, 1.0, 1.0 }; - - vector_solver.CalculateRateConstants(vector_state); - - auto vector_result = vector_solver.Solve(2 * H, vector_state); - EXPECT_NE(vector_result.stats_.singular_, 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 56b65131c..cfe106d09 100644 --- a/test/unit/solver/test_linear_solver_policy.hpp +++ b/test/unit/solver/test_linear_solver_policy.hpp @@ -134,13 +134,12 @@ void testDenseMatrix() auto lu = micm::LuDecomposition::GetLUMatrices(A, 0); auto lower_matrix = std::move(lu.first); auto upper_matrix = std::move(lu.second); - bool is_singular = false; // Only copy the data to the device when it is a CudaMatrix CopyToDeviceSparse(lower_matrix); CopyToDeviceSparse(upper_matrix); - solver.Factor(A, lower_matrix, upper_matrix, is_singular); + solver.Factor(A, lower_matrix, upper_matrix); solver.template Solve(x, lower_matrix, upper_matrix); // Only copy the data to the host when it is a CudaMatrix @@ -188,13 +187,12 @@ void testRandomMatrix(std::size_t number_of_blocks) auto lu = micm::LuDecomposition::GetLUMatrices(A, 0); auto lower_matrix = std::move(lu.first); auto upper_matrix = std::move(lu.second); - bool is_singular = false; // Only copy the data to the device when it is a CudaMatrix CopyToDeviceSparse(lower_matrix); CopyToDeviceSparse(upper_matrix); - solver.Factor(A, lower_matrix, upper_matrix, is_singular); + solver.Factor(A, lower_matrix, upper_matrix); solver.template Solve(x, lower_matrix, upper_matrix); // Only copy the data to the host when it is a CudaMatrix @@ -254,13 +252,12 @@ void testExtremeInitialValue(std::size_t number_of_blocks, double initial_value) auto lu = micm::LuDecomposition::GetLUMatrices(A, initial_value); auto lower_matrix = std::move(lu.first); auto upper_matrix = std::move(lu.second); - bool is_singular = false; // Only copy the data to the device when it is a CudaMatrix CopyToDeviceSparse(lower_matrix); CopyToDeviceSparse(upper_matrix); - solver.Factor(A, lower_matrix, upper_matrix, is_singular); + solver.Factor(A, lower_matrix, upper_matrix); // Only copy the data to the host when it is a CudaMatrix CopyToHostDense(lower_matrix); @@ -308,13 +305,12 @@ void testDiagonalMatrix(std::size_t number_of_blocks) auto lu = micm::LuDecomposition::GetLUMatrices(A, 0); auto lower_matrix = std::move(lu.first); auto upper_matrix = std::move(lu.second); - bool is_singular = false; // Only copy the data to the device when it is a CudaMatrix CopyToDeviceSparse(lower_matrix); CopyToDeviceSparse(upper_matrix); - solver.Factor(A, lower_matrix, upper_matrix, is_singular); + solver.Factor(A, lower_matrix, upper_matrix); solver.template Solve(x, lower_matrix, upper_matrix); // Only copy the data to the host when it is a CudaMatrix diff --git a/test/unit/solver/test_lu_decomposition.cpp b/test/unit/solver/test_lu_decomposition.cpp index 13b089b47..190c46a5b 100644 --- a/test/unit/solver/test_lu_decomposition.cpp +++ b/test/unit/solver/test_lu_decomposition.cpp @@ -18,11 +18,6 @@ TEST(LuDecomposition, DenseMatrixStandardOrdering) testDenseMatrix(); } -TEST(LuDecomposition, SingularMatrixStandardOrdering) -{ - testSingularMatrix(); -} - TEST(LuDecomposition, RandomMatrixStandardOrdering) { testRandomMatrix(1); @@ -51,14 +46,6 @@ TEST(LuDecomposition, DenseMatrixVectorOrdering) testDenseMatrix(); } -TEST(LuDecomposition, SingluarMatrixVectorOrdering) -{ - testSingularMatrix(); - testSingularMatrix(); - testSingularMatrix(); - testSingularMatrix(); -} - TEST(LuDecomposition, RandomMatrixVectorOrdering) { testRandomMatrix(5); diff --git a/test/unit/solver/test_lu_decomposition_policy.hpp b/test/unit/solver/test_lu_decomposition_policy.hpp index 9c9067a00..f1657c168 100644 --- a/test/unit/solver/test_lu_decomposition_policy.hpp +++ b/test/unit/solver/test_lu_decomposition_policy.hpp @@ -122,33 +122,11 @@ void testDenseMatrix() LuDecompositionPolicy lud = LuDecompositionPolicy(A); auto LU = micm::LuDecomposition::GetLUMatrices(A, 0); - bool is_singular{ false }; - lud.template Decompose(A, LU.first, LU.second, is_singular); + lud.template Decompose(A, LU.first, LU.second); check_results( A, LU.first, LU.second, [&](const double a, const double b) -> void { EXPECT_NEAR(a, b, 1.0e-10); }); } -template -void testSingularMatrix() -{ - SparseMatrixPolicy A = SparseMatrixPolicy( - SparseMatrixPolicy::Create(2).InitialValue(0).WithElement(0, 0).WithElement(0, 1).WithElement(1, 0).WithElement(1, 1)); - - A[0][0][0] = 0; - A[0][0][1] = 1; - A[0][1][0] = 1; - A[0][1][1] = 1; - - LuDecompositionPolicy lud = LuDecompositionPolicy(A); - auto LU = micm::LuDecomposition::GetLUMatrices(A, 1.0E-30); - bool is_singular{ false }; - lud.template Decompose(A, LU.first, LU.second, is_singular); - EXPECT_TRUE(is_singular); - A[0][0][0] = 12; - lud.template Decompose(A, LU.first, LU.second, is_singular); - EXPECT_FALSE(is_singular); -} - template void testRandomMatrix(std::size_t number_of_blocks) { @@ -171,8 +149,7 @@ void testRandomMatrix(std::size_t number_of_blocks) LuDecompositionPolicy lud = LuDecompositionPolicy(A); auto LU = micm::LuDecomposition::GetLUMatrices(A, 0); - bool is_singular{ false }; - lud.template Decompose(A, LU.first, LU.second, is_singular); + lud.template Decompose(A, LU.first, LU.second); check_results( A, LU.first, LU.second, [&](const double a, const double b) -> void { EXPECT_NEAR(a, b, 1.0e-9); }); } @@ -213,9 +190,7 @@ void testExtremeValueInitialization(std::size_t number_of_blocks, double initial CopyToDevice(LU.first); CopyToDevice(LU.second); - bool is_singular{ false }; - - lud.template Decompose(A, LU.first, LU.second, is_singular); + lud.template Decompose(A, LU.first, LU.second); CopyToHost(LU.first); CopyToHost(LU.second); @@ -241,8 +216,7 @@ void testDiagonalMatrix(std::size_t number_of_blocks) LuDecompositionPolicy lud = LuDecompositionPolicy(A); auto LU = micm::LuDecomposition::GetLUMatrices(A, 0); - bool is_singular{ false }; - lud.template Decompose(A, LU.first, LU.second, is_singular); + lud.template Decompose(A, LU.first, LU.second); check_results( A, LU.first, LU.second, [&](const double a, const double b) -> void { EXPECT_NEAR(a, b, 1.0e-10); }); } \ No newline at end of file diff --git a/test/unit/solver/test_rosenbrock.cpp b/test/unit/solver/test_rosenbrock.cpp index 06f3cd4bf..357f2b683 100644 --- a/test/unit/solver/test_rosenbrock.cpp +++ b/test/unit/solver/test_rosenbrock.cpp @@ -134,104 +134,4 @@ TEST(RosenbrockSolver, VectorNormalizedError) testNormalizedErrorDiff(VectorBuilder<4>(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()), 3); testNormalizedErrorDiff(VectorBuilder<8>(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()), 5); testNormalizedErrorDiff(VectorBuilder<10>(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()), 3); -} - -TEST(RosenbrockSolver, SingularSystemZeroInBottomRightOfU) -{ - auto params = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters(); - params.check_singularity_ = true; - auto standard = StandardBuilder(params); - auto vector = VectorBuilder<3>(params); - - auto standard_solver = getSingularSystemZeroInBottomRightOfU(standard).SetNumberOfGridCells(1).Build(); - auto vector_solver = getSingularSystemZeroInBottomRightOfU(vector).SetNumberOfGridCells(4).Build(); - - auto standard_state = standard_solver.GetState(); - auto vector_state = vector_solver.GetState(); - - double k1 = -2; - double k2 = 1.0; - - standard_state.SetCustomRateParameter("r1", k1); - standard_state.SetCustomRateParameter("r2", k2); - - vector_state.SetCustomRateParameter("r1", { k1, k1, k1, k1 }); - vector_state.SetCustomRateParameter("r2", { k2, k2, k2, k2 }); - - standard_state.variables_[0] = { 1.0, 1.0 }; - vector_state.variables_[0] = { 1.0, 1.0 }; - vector_state.variables_[1] = { 1.0, 1.0 }; - vector_state.variables_[2] = { 1.0, 1.0 }; - vector_state.variables_[3] = { 1.0, 1.0 }; - - // to get a jacobian with an LU factorization that contains a zero on the diagonal - // of U, we need det(alpha * I - jacobian) = 0 - // for the system above, that means we have to have alpha + k1 + k2 = 0 - // in this case, one of the reaction rates will be negative but it's good enough to - // test the singularity check - // alpha is 1 / (H * gamma), where H is the time step and gamma is the gamma value from - // the rosenbrock paramters - // so H needs to be 1 / ( (-k1 - k2) * gamma) - // since H is positive we need -k1 -k2 to be positive, hence the smaller, negative value for k1 - double H = 1 / ((-k1 - k2) * params.gamma_[0]); - standard_solver.solver_.parameters_.h_start_ = H; - vector_solver.solver_.parameters_.h_start_ = H; - - standard_solver.CalculateRateConstants(standard_state); - vector_solver.CalculateRateConstants(vector_state); - - auto standard_result = standard_solver.Solve(2 * H, standard_state); - EXPECT_NE(standard_result.stats_.singular_, 0); - - auto vector_result = vector_solver.Solve(2 * H, vector_state); - EXPECT_NE(vector_result.stats_.singular_, 0); -} - -TEST(RosenbrockSolver, SingularSystemZeroAlongDiagonalNotBottomRight) -{ - auto params = micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters(); - - double k1 = -1.0; - double k2 = -1.0; - double k3 = 1.0; - - // to get a jacobian with an LU factorization that contains a zero on the diagonal - // of U, we need det(alpha * I - jacobian) = 0 - // for the system above, that means we have to set alpha = -k1, or alpha=-k2, or alpha=k3 - double H = 1 / (-k1 * params.gamma_[0]); - - params.check_singularity_ = true; - params.h_start_ = H; - - auto standard = StandardBuilder(params); - auto vector = VectorBuilder<3>(params); - - auto standard_solver = getSolverForSingularSystemOnDiagonal(standard).SetNumberOfGridCells(1).Build(); - auto vector_solver = getSolverForSingularSystemOnDiagonal(vector).SetNumberOfGridCells(4).Build(); - - auto standard_state = standard_solver.GetState(); - auto vector_state = vector_solver.GetState(); - - standard_state.SetCustomRateParameter("r1", k1); - standard_state.SetCustomRateParameter("r2", k2); - standard_state.SetCustomRateParameter("r3", k3); - - vector_state.SetCustomRateParameter("r1", { k1, k1, k1, k1 }); - vector_state.SetCustomRateParameter("r2", { k2, k2, k2, k2 }); - vector_state.SetCustomRateParameter("r3", { k3, k3, k3, k3 }); - - standard_state.variables_[0] = { 1.0, 1.0, 1.0 }; - vector_state.variables_[0] = { 1.0, 1.0, 1.0 }; - vector_state.variables_[1] = { 1.0, 1.0, 1.0 }; - vector_state.variables_[2] = { 1.0, 1.0, 1.0 }; - vector_state.variables_[3] = { 1.0, 1.0, 1.0 }; - - standard_solver.CalculateRateConstants(standard_state); - vector_solver.CalculateRateConstants(vector_state); - - auto standard_result = standard_solver.Solve(2 * H, standard_state); - EXPECT_NE(standard_result.stats_.singular_, 0); - - auto vector_result = vector_solver.Solve(2 * H, vector_state); - EXPECT_NE(vector_result.stats_.singular_, 0); -} +} \ No newline at end of file diff --git a/test/unit/solver/test_rosenbrock_solver_policy.hpp b/test/unit/solver/test_rosenbrock_solver_policy.hpp index 855d91800..daf37714a 100644 --- a/test/unit/solver/test_rosenbrock_solver_policy.hpp +++ b/test/unit/solver/test_rosenbrock_solver_policy.hpp @@ -96,68 +96,4 @@ void testAlphaMinusJacobian(SolverBuilderPolicy builder, std::size_t number_of_g EXPECT_NEAR(jacobian[i_cell][4][2], -53.6, 1.0e-5); EXPECT_NEAR(jacobian[i_cell][4][4], 42.042 - 1.0, 1.0e-5); } -} - -template -SolverBuilderPolicy getSingularSystemZeroInBottomRightOfU(SolverBuilderPolicy builder) -{ - // A -> B - // B -> A - - auto a = micm::Species("a"); - auto b = micm::Species("b"); - - micm::Phase gas_phase{ std::vector{ a, b } }; - - micm::Process r1 = micm::Process::Create() - .SetReactants({ a }) - .SetProducts({ Yields(b, 1) }) - .SetPhase(gas_phase) - .SetRateConstant(micm::UserDefinedRateConstant({ .label_ = "r1" })); - - micm::Process r2 = micm::Process::Create() - .SetReactants({ b }) - .SetProducts({ Yields(a, 1) }) - .SetPhase(gas_phase) - .SetRateConstant(micm::UserDefinedRateConstant({ .label_ = "r2" })); - - return builder.SetSystem(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase })) - .SetReactions(std::vector{ r1, r2 }) - .SetReorderState(false); -} - -template -SolverBuilderPolicy getSolverForSingularSystemOnDiagonal(SolverBuilderPolicy builder) -{ - // A -> B, k1 - // B -> C, k2 - // C -> A, k3 - - auto a = micm::Species("a"); - auto b = micm::Species("b"); - auto c = micm::Species("c"); - - micm::Phase gas_phase{ std::vector{ a, b, c } }; - - micm::Process r1 = micm::Process::Create() - .SetReactants({ a }) - .SetProducts({ Yields(b, 1) }) - .SetPhase(gas_phase) - .SetRateConstant(micm::UserDefinedRateConstant({ .label_ = "r1" })); - - micm::Process r2 = micm::Process::Create() - .SetReactants({ b }) - .SetProducts({ Yields(c, 1) }) - .SetPhase(gas_phase) - .SetRateConstant(micm::UserDefinedRateConstant({ .label_ = "r2" })); - - micm::Process r3 = micm::Process::Create() - .SetReactants({ c }) - .SetProducts({ Yields(a, 1) }) - .SetPhase(gas_phase) - .SetRateConstant(micm::UserDefinedRateConstant({ .label_ = "r3" })); - - return builder.SetSystem(micm::System(micm::SystemParameters{ .gas_phase_ = gas_phase })) - .SetReactions(std::vector{ r1, r2, r3 }) - .SetReorderState(false); } \ No newline at end of file