SUNDIALS: update ARKode to version 5.4.0 #11278
Conversation
There was a problem hiding this comment.
Thanks for the contribution. I will go through the code later today (to learn something regarding Sundials) and then let the actual review do by someone who is actually using Sundials!
In the meantime: maybe you could associate your LNM-email address with your GitHub-account. Then GitHub would be able to pick up the information that you have made the commits ;)
There was a problem hiding this comment.
I left some comments. Overall it looks already very good! Thanks a lot for the affords, @sebproell!
There are some issues which need to be discussed/addressed by all of us:
- There are many TODOs regarding different options. @sebproell: Do you plan to provide more options? I am fine if we have for now only a limited number options. However, we should document the limitations and maybe add at some places some asserts if needed.
- What is the problem regarding the vector? Does ARKode want its own vector type? Could you elaborate the issues and what you mean by "this method really doesn't belong here and should be removed once N_Vector "understands" our vectors"? Do you have any idea how to make it work without the copy step? How could we remove the
create_vector()function that takes a dummy vector? Furthermore, could we cache the vectors created bycreate_vector: it seem like a waste of operation to create a new vector eachvmult()? - IDA and KINSOL are not implemented currently for newest version. Till the next release we should make these also work. @sebproell Any ideas how much has changed there!? Maybe there is volunteer who would like to take up the task in a follow-up PR?
source/sundials/arkode.cc
Outdated
| // TODO currently no distinction between left and right preconditioning is | ||
| // possible | ||
| int status = p_solve_fn(P_data, sun_src, sun_dst, tol, 0); |
There was a problem hiding this comment.
Why not? Where does:
dealii/include/deal.II/sundials/arkode.h
Lines 1068 to 1069 in 6a03dd0
come into play?
There was a problem hiding this comment.
I couldn't find anything on left and right preconditioning in the dealii docs. Therefore, I decided to not use this distinction if the user provides a custom solver. The line you are commenting on only becomes relevant when the preconditioner is passed back to a user-written solver function. I added some more remarks in the doc string.
|
Let me make some remarks regarding SUNDIALS' They define N_Vector as typedef struct _generic_N_Vector *N_Vector;
struct _generic_N_Vector {
void *content;
struct _generic_N_Vector_Ops *ops;
};Users are supposed to attach their actual vector to the I think in the future this behavior should be replaced with an actual N_Vector wrapper for the supported dealii-vectors. I can do this in a follow-up PR, in fact I already started a first try in this direction. If this is to be changed anyway I am not sure how much effort should be spent on the current implementation. The ugly @peterrum I think this also answers a couple of your inline remarks (thanks btw). |
|
First thank you @sebproell for tackling such a big topic in the first PR. @sebproell and I have discussed some of the design choices offline, so I will later mainly go through the details and invite others to comment on the new interface as I am biased. I think one of the issues to decide collectively regarding this PR is the question of compatibility: Do we want to support both the old version of SUNDIALS (3.2.1 IIRC) and the new one provided with wrappers by this PR, or do we keep both versions? This PR nicely packages things into |
include/deal.II/sundials/arkode.h
Outdated
| @@ -49,6 +52,24 @@ | |||
|
|
|||
| DEAL_II_NAMESPACE_OPEN | |||
|
|
|||
| // Forward declarations | |||
| # if DEAL_II_SUNDIALS_VERSION_GTE(5, 4, 0) | |||
There was a problem hiding this comment.
What do we do for versions between 4 and 5.4? Are these all classes that genuinely only showed up with 5.4, or have they been around since 4.0?
There was a problem hiding this comment.
Sam of course in similar places below
There was a problem hiding this comment.
Actually the linear and nonlinear solver interfaces were introduced with 4.0.0. I adapted the implementation (will push soon) to also compile with SUNDIALS 4.0.0 and had to adapt some tests. Unfortunately, one convenient function (SUNLinSolNewEmpty) is not yet available there but it is basically necessary.
I only had a brief look. It looks like we also relied on some things in an implementation header that is no longer installed with SUNDIALS. Again the changes should mainly concern the linear and nonlinear solvers. I would assume a similar strategy as for ARKODE might be applicable. |
| * This file contains a verbatim copy of code snippets distributed within the | ||
| * SUNDIALS package, see the license below. | ||
| */ |
There was a problem hiding this comment.
As already mentioned I had to copy two functions from the newest SUNDIALS version to use 4.0.0. I put them in this file since I have no idea how to handle such a case especially with regards to licenses.
Could you point out the correct place in the documentation that specifies which operations have to be defined. Furthermore, if you already have started to implement a wrapper, may I ask you to share it (either by pointing out the right branch or by already creating a WIP PR). Thx. |
See https://computing.llnl.gov/sites/default/files/public/ark_guide-4.5.0_0.pdf starting on page 307. And yes, I will share a first shot at implementing a few operations later. |
See sebproell#1. I implemented some basic clone and destroy operations and the linear sum. |
There was a problem hiding this comment.
I have a few cosmetic changes but it looks great otherwise. Given that this is @sebproell's first PR, I think we can leave some of the cleanup and extension to other SUNDIALS packages besides ARKode to later PRs. @luca-heltai given that you have been involved in the original wrapper via #5156, would you mind having a look? I suggest we try to get this PR ready in the course of the next days, given that @sebproell has some follow-up work pending.
|
/rebuild |
|
@sebproell could you please squash the commits into one or a few logically separate commits and then force-push the branch again? |
sebproell
force-pushed
the
update-sundials-arkode-5.4.0
branch
from
062781c
to
06cbef4
Compare
Provide new std::functions to interface with SUNDIALS arkode module IDA and KINSOL not updated
sebproell
force-pushed
the
update-sundials-arkode-5.4.0
branch
from
06cbef4
to
7ec3423
Compare
|
I will try to review this patch in the next two or three days. In the mean time: thanks @sebproell for your work! :) |
There was a problem hiding this comment.
I really appreciate the amount of work you put into this, however I see a few (big) issues with this PR that we need to address (in followup PR) but BEFORE 9.3:
- it keeps intact the old interface to SUNDIALS 3.2, without really exploiting the new structure of SUNDIALS > 5 and of deal.II 9.2
- it exposes in several ways the internals of SUNDIALS, and it adds essentially one-to-one translations to C code in ARKODE, without exploiting the new concepts
- it replicates concepts we already have, using SUNDIALS notations
The new SUNDIALS interface has gone through a lot of work to separate linear solvers, non linear solvers, and linear operators. In particular now all of sundials interfaces use the concepts of
- NVector
- SUNMatrix
- SUNLinearSolver
- SUNNonlinearSolver
The version that deal.II supported, only had a uniform use of the NVector. Every package (ARKODE, IDA, KINSOL) had slightly different solvers, internal structures, nonlinear solvers, etc.
Introducing support to SUNDIALS 5.0 without interfaces to
- SUNMatrix
- SUNLinearSolver
- SUNNonlinearSolver
Will make our life very difficult in the future (and it already made it difficult for you: now you had to workaround the fact that providing a custom solver is different in this new version). The undocumented way to supply a solver is now supported (and documented): it is called SUNLinearSolver.
Providing a wrapper for ARKODE (or IDA, or CVODE) without providing a wrapper for SUNLinearSolver and SUNNonlinearSolver basically produces a one-to-one copy of SUNDIALS objects inside our library, without really simplifying our usage of SUNDIALS.
The original design of IDA and ARKODE interface in deal.II was that any deal.II user that is doing serious time dependent problem has a mass matrix (free) at his disposal, a matrix assembly (or matrix vector product evaluation) routine, and most likely a solver (including a preconditioner), since PDE problems are difficult to assemble, and difficult to solve without providing good preconditioners and solvers.
In other words, we almost always have at our disposal a concept of SUNLinearSolver, a SUNMatrix and possibly several NVector already.
The approach of the original interface we had was to plug our solvers into SUNDIALS (and not to expose sundials solvers to us: PETSc, Trilinos, and deal.II itself do a very good job with solvers/preconditioners for our use case. SUNDIALS ones were not designed with PDEs in mind, and there was no need to expose them inside deal.II).
This is now possible by design: if we provide an interface to SUNLinSolver, then IDA/ARKODE/CVODE, etc. can use that.
I suspect that adapting the old design to the new SUNDIALS version without exploiting this is a very bad decision.
Having said this: you already did 95% of the job. But this will likely need to be restructured (a lot!) to really exploit the new SUNDIALS interface. I suggest we move forward and merge this PR. Next step should be to coordinate together (I had started working on the interface already, but then teaching started, and with it, my time with deal.II dropped) and make sure we extract common parts into a separate header file (like the create_nvector function), and provide interfaces to our own LinearOperator + Solvers as SUNLinSolver and SUNLinOperator, so that we can reuse all of them in all of SUNDIALS packages.
include/deal.II/sundials/arkode.h
Outdated
| @@ -907,6 +1305,62 @@ namespace SUNDIALS | |||
| # endif // DEAL_II_WITH_PETSC | |||
| }; | |||
|
|
|||
| # if DEAL_II_SUNDIALS_VERSION_GTE(5, 4, 0) | |||
There was a problem hiding this comment.
This structure is identical to a dealii LinearOperator (with very small differences). Do you think it would be possible to refactor our interface to arkode to actually take two LinearOperator objects?
There was a problem hiding this comment.
Yes, I read through the docs of LinearOperator and realized that the vmult is a public field that can just be set externally. Then my classes are unnecessary. Also should answer your question below.
| * built internally by SUNDIALS based on user settings. The parameters are | ||
| * interpreted as follows: | ||
| * | ||
| * @param[in] op A LinearOperator that applies the matrix vector product |
There was a problem hiding this comment.
The documentation confuses me. As far as I can understand, you introduced two new types of LinearOperator objects (SundialsOperator, and SundialsPreconditioner) which are essentially copies of our own LinearOperator type. The documentation here states what should have been the right interface (i.e.: use two of our dealii::LinearOperator objects), but then the signature is different, and requires two new objects.
Notice that the original interfaces to sundials where constructed before the introduction of the LinearOperator concept. Now that we have that concept, and that Sundials moved in the same direction, it should in principle be sufficient to pass two linear operators.
The general strategy, with the SUNDIALS wrappers, was centered around the fact that, as deal.II users and developers, we almost always build our own matrices, preconditioners, and solvers. The interface to SUNDIALS was supposed to be "a few lambda away" from existing codes.
Wouldn't it have been better to allow sundials to interpret our LinearOperator classes as objects that can be used directly inside sundials, instead of introducing (yet another two) interpretation of the same object?
| @@ -355,6 +435,8 @@ namespace SUNDIALS | |||
| const unsigned int maximum_non_linear_iterations = 10, | |||
| const bool implicit_function_is_linear = false, | |||
| const bool implicit_function_is_time_independent = false, | |||
| const bool mass_is_time_independent = false, | |||
| const int anderson_acceleration_subspace = 3, | |||
There was a problem hiding this comment.
| const int anderson_acceleration_subspace = 3, | |
| const int anderson_acceleration_subspace_size = 3, |
| @@ -370,41 +452,15 @@ namespace SUNDIALS | |||
| , implicit_function_is_linear(implicit_function_is_linear) | |||
| , implicit_function_is_time_independent( | |||
| implicit_function_is_time_independent) | |||
| , mass_is_time_independent(mass_is_time_independent) | |||
| , anderson_acceleration_subspace(anderson_acceleration_subspace) | |||
There was a problem hiding this comment.
| , anderson_acceleration_subspace(anderson_acceleration_subspace) | |
| , anderson_acceleration_subspace_size(anderson_acceleration_subspace_size) |
| @@ -431,6 +487,9 @@ namespace SUNDIALS | |||
| implicit_function_is_linear); | |||
| prm.add_parameter("Implicit function is time independent", | |||
| implicit_function_is_time_independent); | |||
| prm.add_parameter("Mass is time independent", mass_is_time_independent); | |||
| prm.add_parameter("Anderson-acceleration subspace", | |||
There was a problem hiding this comment.
| prm.add_parameter("Anderson-acceleration subspace", | |
| prm.add_parameter("Anderson-acceleration subspace size", |
| @@ -431,6 +487,9 @@ namespace SUNDIALS | |||
| implicit_function_is_linear); | |||
| prm.add_parameter("Implicit function is time independent", | |||
| implicit_function_is_time_independent); | |||
| prm.add_parameter("Mass is time independent", mass_is_time_independent); | |||
| prm.add_parameter("Anderson-acceleration subspace", | |||
| anderson_acceleration_subspace); | |||
There was a problem hiding this comment.
| anderson_acceleration_subspace); | |
| anderson_acceleration_subspace_size); |
| * Number of subspace vectors to use for Anderson acceleration. Only | ||
| * meaningful if the packaged SUNDIALS fixed-point solver is used. | ||
| */ | ||
| int anderson_acceleration_subspace; |
There was a problem hiding this comment.
| int anderson_acceleration_subspace; | |
| int anderson_acceleration_subspace_size; |
| @@ -541,6 +612,44 @@ namespace SUNDIALS | |||
| void | |||
| reset(const double t, const double h, const VectorType &y); | |||
|
|
|||
| /*! | |||
There was a problem hiding this comment.
We had long discussions about exposing SUNDIALS internal vector, solvers, and matrix types the first time around when we introduced the SUNDIALS wrappers, and the general consensus was that none of this should have been exposed to users. The original approach was to copy back and forth between vector types, with the agreement that whenever possible, we should have transformed those copies into views. We never got around to actually performing those copies into view only actions, partially because of the different layouts of the MPI versions of the vectors, and because there was no easy way to perform the conversion on the other way around (i.e., deal.II does not offer yet a View concept into its own vector types (even thought it used to, with the class VectorView, only in serial).
Since you are still performing these copies internally, why do you need to expose this interface externally?
We have a deal.II/sundials/copy.h file (that could be renamed to vector.h), that provides all these functions already. If you prefer/need this interface, you should add it directly in that header file, and remove it from arkode.h. Otherwise we'd have to do the same for IDA/CVODE, etc.
There was a problem hiding this comment.
I proposed some changes in the "view instead of copy" direction in #11395.
| * functions | ||
| */ | ||
| void * | ||
| get_arkode_memory() const; |
There was a problem hiding this comment.
Please return an ARKodeMem, and before returning it, Assert that the conversion is successfull. Or do you have a reason to keep it void *? If we know it is ARKodeMem, let's return an ARKodeMem type.
|
@luca-heltai Thanks for your input and feedback! Overall, I agree with all your points and @peterrum already told me I should have written a little more of this in the introduction of this PR. Will do next time. I am planning to use SUNDIALS in my projects and wanted a quick update to the newest version. I am happy to help in further efforts to update and unify the different modules. see #11395 for some first steps Should I address your remarks about the current implementation here or in a followup? |
|
I'll merge this, and open an issue to track of the left overs. Thanks for all your work and your patience! |
This PR proposes a new set of std::functions to hook into a newer version of SUNDIALS ARKode package. Quite some changes happened there compared to the last supported version 3.2.1. In particular, the old way of specifying linear solvers (via some undocumented internal headers) is no longer possible.
Since this is my first PR here I am counting on @kronbichler and @peterrum to mention some appropriate reviewers or add necessary information.