-
Notifications
You must be signed in to change notification settings - Fork 3
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Vector-ordered Rosenbrock Solver #172
Conversation
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
looks good! just a couple questions and suggestions. the only important one is that that the number of grid cells won't be known at compile time
if(MUSICA_ENABLE_OPENMP) | ||
target_link_libraries(test_${TEST_NAME} PUBLIC OpenMP::OpenMP_CXX OpenMP::OpenMP_Fortran) | ||
target_link_libraries(test_${TEST_NAME} PUBLIC OpenMP::OpenMP_CXX) | ||
endif() |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Hi @mattldawson and @K20shores, what do you think about this change? Matt and I briefly discussed about this. The issue is that OpenMP::OpenMP_Fortran
is not found when test_micm_c_api
is built. I removed it because this line belongs to function(create_standard_test_cxx)
and I assumed that it shouldn't need OpenMp_Fortran
for c++ tests.
Matt mentioned we could specify fortran as one of the languages in the top-level CMakeLists.txt
and update find_package
in cmake/dependencies.cmake
project(musica-distribution VERSION 0.7.0 LANGUAGES C CXX Fortran)
find_package(OpenMP REQUIRED COMPONENTS C CXX Fortran)
I think that makes sense although I think it is useful to separate the fortran stuff from MUSICA-C
as it is since we eventually pursue fortran-free. I'd like to hear what you all think.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I would like to keep fortran separated, but with the inclusion of tuvx into the musica-c library that line is blurred. Because of that I think we should go with Matt's suggestion. The other option is to check if TUVX is enabled before we check for OpenMP. If it is, enable the fortran language so that the fortran part of OpenMP is found, if that's possible
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I created a separate issue to tackle this.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think it all looks fine as is, but we could make it slightly better. I wonder if we could remove the need to have two member variable pointers to different solver types and two different build functions when only the matrix types are different
! We could use Fortran 2023 enum type feature if Fortran 2023 is supported | ||
! https://fortran-lang.discourse.group/t/enumerator-type-in-bind-c-derived-type-best-practice/5947/2 | ||
enum, bind(c) | ||
enumerator :: Rosenbrock = 1 | ||
enumerator :: RosenbrockStandardOrder = 2 | ||
end enum |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
As long as our CI tests pass then I'm fine with this. The only potential issue maybe is the NAG compiler and if it supports this and if we need to compile with NAG ever.
But I like this as it is
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Jesse mentioned that it's probably not yet to use Fortran 2023 features but I'll leave the comment so we can revisit later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I like this too
/// @brief Vector-ordered Rosenbrock solver type | ||
using DenseMatrixVector = micm::VectorMatrix<double, MICM_VECTOR_MATRIX_SIZE>; | ||
using SparseMatrixVector = micm::SparseMatrix<double, micm::SparseMatrixVectorOrdering<MICM_VECTOR_MATRIX_SIZE>>; | ||
using RosenbrockVectorType = typename micm::RosenbrockSolverParameters:: | ||
template SolverType<micm::ProcessSet, micm::LinearSolver<SparseMatrixVector, micm::LuDecomposition>>; | ||
using Rosenbrock = micm::Solver<RosenbrockVectorType, micm::State<DenseMatrixVector, SparseMatrixVector>>; | ||
std::unique_ptr<Rosenbrock> rosenbrock_; | ||
|
||
std::unique_ptr<Rosenbrock> solver_; | ||
/// @brief Standard-ordered Rosenbrock solver type | ||
using DenseMatrixStandard = micm::Matrix<double>; | ||
using SparseMatrixStandard = micm::SparseMatrix<double, micm::SparseMatrixStandardOrdering>; | ||
using RosenbrockStandardType = typename micm::RosenbrockSolverParameters:: | ||
template SolverType<micm::ProcessSet, micm::LinearSolver<SparseMatrixStandard, micm::LuDecomposition>>; | ||
using RosenbrockStandard = micm::Solver<RosenbrockStandardType, micm::State<DenseMatrixStandard, SparseMatrixStandard>>; | ||
std::unique_ptr<RosenbrockStandard> rosenbrock_standard_; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Are we allowed to have a pointer to a micm::Solver
without specifying the template types?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
or maybe we use std::variant
?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Do you mean something like std::unique_ptr<micm::Solver> solver
? If it is the case, I tried it but it didn't work, (template argument is invalid
). UnlessMICM
is a template class which is not a viable option, I think all the types must be specified.
I tried implementing std::variant
and realized that you still need to do type check to get the data that variant holds, and the syntax can be quite tricky. For example, MICM::Solve( lots of params ... )
// Current
micm::State state = solver->GetState();
// std::variant
if(std::get_if<Rosenbrock>(&solver))
micm::State state = std::get<Rosenbrock>(&solver)->GetState();
if(std::get_if<RosenbrockStandardOrder>(&solver))
micm::State state = std::get<RosenbrockStandardOrder>(&solver)->GetState();
What do you think? The fact that C doesn't support any form of overloading, we quite have limitations here.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Hmm. Good point. What you have seems fine anyways
micm::SolverBuilder< | ||
micm::RosenbrockSolverParameters, | ||
micm::VectorMatrix<double, MICM_VECTOR_MATRIX_SIZE>, | ||
micm::SparseMatrix<double, micm::SparseMatrixVectorOrdering<MICM_VECTOR_MATRIX_SIZE>>, | ||
micm::ProcessSet, | ||
micm::LinearSolver< | ||
micm::SparseMatrix<double, micm::SparseMatrixVectorOrdering<MICM_VECTOR_MATRIX_SIZE>>, | ||
micm::LuDecomposition>>(micm::RosenbrockSolverParameters::ThreeStageRosenbrockParameters()) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
maybe we should make something similar to micm::CpuSolverBuilder
but make it be micm::VectorizedCpuSolverBuilder
which is templated on the on vector size. Then it would be able to do all of this gross stuff for us. As an added benefit we should be able to add some alias types inside all the solver builder aliases that give us the type of the solver so that inside the header file we would only need to do something like
using builderType.= micm::CpuSolverBuilder;
using solverType = builderType::RosenbrockType;
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I like this idea. I created an issue for this in MICM. I think this can be a separate issue because alias-declaration for builder type should be done first at MICM.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
looks great!
Closes #55
Update
musica
build optionsDoesn't support multi grid cells and vector matrix size greater than 1. (the way assigning the concentrations in Solve function)
I decided to submit this with the limitation described above because this task currently blocks implementing Backward Euler as well as a separate issue was created to support multiple grid cells.
These comments are not valid anymore but left for record to describe where the bug found.
Added compilation flag for the number of grid cells:
MUSICA_SET_NUM_GRID_CELLS
For vector-ordered matrix solver, when number of grid cell is less than
L
value, for example,num_grid_cell=1
,L=4
,num_species=5
, it fails during solving because the Copy() call throws an error indicating that the dimension is not the same; copying other vector to this vector(5 /= 20)
. SinceL
is decided at compile time, checks should be included in the build script