Introduce SUNDIALS N_Vector module #11395

sebproell · 2020-12-18T14:35:10Z

This PR proposes a dealii wrapper for the SUNDIALS' N_Vector module as part of #11380. All of the functionality is in internal because the conversion between vectors is not of interest for a user of a SUNDIALS integrator. I tried to further separate the only relevant functions for an author of the ARKODE, IDA, etc. wrapper into n_vector.h. The functions are: NVectorView, a class to view a VectorType as an N_Vector and unwrap_nvector() and unwrap_nvector_const() which are used to retrieve the internally stored dealii vector from an N_Vector.

This PR doesn't yet make use of the new functionality in any time integrator. I will wait for #11278 to do so but I have already checked that it works (for ARKode at least). Also, there are tests for all the vector operations that are already supported.

Commit b3f1761 contains the relevant changes. The later commits are to (hopefully) fix the tests. I suggest to rebase when #11278 is merged. Note that the tests are not actually run here since they currently require SUNDIALS > 5.0.0. I can later backport functionality or, if support for older versions will be ended in the foreseeable future, we can keep the old copying mechanism for those versions.

sebproell

I highlighted some comments and questions

sebproell · 2020-12-18T14:37:54Z

include/deal.II/sundials/n_vector.h

+     */
+    template <typename VectorType>
+    const VectorType *
+    unwrap_nvector_const(N_Vector v);


Not sure if the calls should explicitly mention the const-ness. Alternative would be to enable something like unwrap_nvector<const VectorType>(...) but this feels unintuitive to me.

sebproell · 2020-12-18T14:42:24Z

include/deal.II/sundials/n_vector.templates.h

+  AssertDimension(x_dealii->size(), z_dealii->size());
+  AssertDimension(x_dealii->size(), y_dealii->size());
+
+  std::transform(x_dealii->begin(),


This is an example of a vector operation that is not easily translated to the dealii vector interface. I went for this naive implementation for now to have a working version in place.

sebproell · 2020-12-18T14:44:36Z

include/deal.II/sundials/n_vector.templates.h

+  //  v->ops->nvwrmsnormmask = undef;
+  //  v->ops->nvmin          = undef;
+  //  v->ops->nvwl2norm      = undef;
+  //  v->ops->nvl1norm       = undef;
+  //  v->ops->nvcompare      = undef;
+  //  v->ops->nvinvtest      = undef;
+  //  v->ops->nvconstrmask   = undef;
+  //  v->ops->nvminquotient  = undef;


For reference, I put some further operations that could be implemented here. In a quick test with the ARKode tests they didn't seem to be necessary.

sebproell · 2020-12-18T14:46:32Z

source/sundials/n_vector.inst.in

+// ---------------------------------------------------------------------
+
+// serial Vector and BlockVector
+for (S : REAL_SCALARS; V : DEAL_II_VEC_TEMPLATES)


I went for a suite of vectors that I thought are representative and use myself. Extension to more vectors should be straight-forward.

sebproell · 2020-12-18T14:52:22Z

tests/sundials/n_vector.cc

+  NVectorView<VectorType> n_vector(vector);
+  auto                    id = N_VGetVectorID(n_vector);


These lines demonstrate the usage of NVectorView. After creating the view, it can be passed to a SUNDIALS function where it implicitly converts to N_Vector. The view is cleaned up when going out of scope.

sebproell · 2020-12-21T11:42:36Z

@kronbichler Can you have a look? This PR is all about dealii vectors, their operations and memory management. Virtually no SUNDIALS knowledge required;)

include/deal.II/sundials/n_vector.templates.h

peterrum

Looks very mature already 🥇 It took some time to understand why so many indirections are needed. One creates a view on deal.II vector NVectorView, which internally creates a NVectorContent. The latter classes can have two states: either it only stores only a pointer to the original vector or it "stores" the actual vector (created by the clone function - or does it?).

Maybe you could provide a small code snippet which shows how unwrap_nvector(), unwrap_nvector_const(), and NVectorView will be used in the ARKODE classes. That might help the next reviewer.

tests/sundials/n_vector.cc

include/deal.II/sundials/n_vector.h

peterrum · 2021-01-10T23:08:07Z

include/deal.II/sundials/n_vector.h

+     * @note N_VDestroy() should not be called on the resulting N_Vector since
+     *   this would lead to a double delete. Let the destructor do this work.


Merge this into the paragraph before. Is there a check?

No and I couldn't come up with an easy one. The problem is that destroy also calls free on the N_Vector.

include/deal.II/sundials/n_vector.h

include/deal.II/sundials/n_vector.templates.h

sebproell

Maybe you could provide a small code snippet which shows how unwrap_nvector(), unwrap_nvector_const(), and NVectorView will be used in the ARKODE classes. That might help the next reviewer.

@peterrum I hope the first three comments in this review help.

In addition, I highlighted a new function make_nvector_view that could be used as a nicer(?) way to construct the view objects. However, it raises some other design questions for the NVectorView

sebproell · 2021-01-12T10:22:43Z

tests/sundials/n_vector.cc

+  auto       vector         = create_test_vector<VectorType>();
+  const auto const_vector   = create_test_vector<VectorType>();


Usage of the new functionality (1/3)
Let's say we have two deal.II vectors vector and const_vector.

sebproell · 2021-01-12T10:23:24Z

tests/sundials/n_vector.cc

+  auto       n_vector       = make_nvector_view(vector);
+  auto       const_n_vector = make_nvector_view(const_vector);


Usage of the new functionality (2/3)
Here, we create two view objects that can be passed to any SUNDIALS function that takes an N_Vector. One would use this approach e.g. to pass an initial solution vector to SUNDIALS.

sebproell · 2021-01-12T10:28:27Z

tests/sundials/n_vector.cc

+  Assert((const_n_vector)->content != nullptr, NVectorTestError());
+
+  // unwrap non-const as non-const
+  auto *vector_unwrapped = unwrap_nvector<VectorType>(n_vector);


Usage of the new functionality (3/3)
This is the oppposite operation: a N_Vector is unwrapped to extract the underlying deal.ii vector. One would use this appraoch e.g. to pass on a dst vector from SUNDIALS to a deal.ii linear solver. Note that after the linear solve no copies or updates must be made since both n_vector and vector_unwrapped refer to the same memory.

sebproell · 2021-01-12T10:37:01Z

include/deal.II/sundials/n_vector.h

+     */
+    template <typename VectorType>
+    NVectorView<VectorType>
+    make_nvector_view(VectorType &vector);


I introduced this new function to utilize template parameter deduction.

sebproell · 2021-01-12T10:39:34Z

include/deal.II/sundials/n_vector.h

+       *     auto view = make_nvector_view(vector);
+       *   @endcode
+       */
+      NVectorView() = default;


See the doc string. I am not 100% sure if this is required. With gcc 7.3.1 it also worked without a default constructor.

Note to myself: seems that this is not required. I will remove it.

sebproell · 2021-01-12T10:41:04Z

include/deal.II/sundials/n_vector.h

+      /**
+       * Constructor. Create view of @p vector.
+       */
+      NVectorView(VectorType &vector);


Could be made private if make_nvector_view is declared as friend. Currently, one could still use the constructor.

sebproell · 2021-01-12T11:26:59Z

It took some time to understand why so many indirections are needed. One creates a view on deal.II vector NVectorView, which internally creates a NVectorContent. The latter classes can have two states: either it only stores only a pointer to the original vector or it "stores" the actual vector (created by the clone function - or does it?).

Correct, to put it differently: the class NVectorContent is the type for objects that are attached to the void* content of the N_Vector type. It manages the memory of the underlying vector if necessary. SUNDIALS creates new vector with N_VClone on an existing vector, which will create an owning NVectorContent. The second class NVectorView is for convenience and takes over creation of a non-owning NVectorContent and attaches it to an empty N_Vector.

Only NVectorView would ever be used in SUNDIALS code by a deal author. A user of deal would ofc never see any of the stuff going on in this PR.

sebproell · 2021-01-21T16:15:11Z

@peterrum @luca-heltai @kronbichler How should we proceed here?

luca-heltai · 2021-01-22T08:25:06Z

/rebuild

luca-heltai · 2021-01-22T08:26:18Z

I think we should merge this, and adjust later if necessary. With this, it should be much easier to interface with sundials, and I think it is a big step forward. I'm in favour of merging early, testing a bit the new interface, and then adjust later if necessary.

sebproell · 2021-01-22T09:59:24Z

Okay! I already have some adaptions to the ARKode integrator to use the new N_Vector module. Will follow up with it soon.

sebproell commented Dec 18, 2020

View reviewed changes

sebproell marked this pull request as ready for review December 18, 2020 16:15

peterrum self-requested a review December 27, 2020 08:03

sebproell mentioned this pull request Dec 29, 2020

SUNDIALS: update ARKode to version 5.4.0 #11278

Merged

luca-heltai mentioned this pull request Dec 29, 2020

SUNDIALS interface. #11435

Open

5 tasks

luca-heltai reviewed Dec 29, 2020

View reviewed changes

include/deal.II/sundials/n_vector.templates.h Show resolved Hide resolved

luca-heltai added this to In Progress in SUNDIALS Interfaces Dec 29, 2020

luca-heltai added this to the Release 9.3 milestone Dec 29, 2020

Introduce SUNDIALS N_Vector module

c527760

sebproell force-pushed the sundials-n_vector branch from ef3cc23 to c527760 Compare January 7, 2021 08:55

peterrum reviewed Jan 11, 2021

View reviewed changes

make_nvector_view and cleanup

eaf7faf

sebproell commented Jan 12, 2021

View reviewed changes

luca-heltai added ready to test Reviewed and ready to merge labels Jan 22, 2021

luca-heltai approved these changes Jan 22, 2021

View reviewed changes

peterrum approved these changes Jan 22, 2021

View reviewed changes

peterrum merged commit c68cbde into dealii:master Jan 22, 2021

sebproell deleted the sundials-n_vector branch January 22, 2021 10:57

sebproell mentioned this pull request Jan 29, 2021

Make N_Vector module compatible with SUNDIALS 3.2.1 #11654

Merged

luca-heltai moved this from In Progress to Done in SUNDIALS Interfaces Jan 30, 2021

This was referenced Jan 30, 2021

SUNDIALS: Use new N_Vector module in ARKode wrapper #11655

Merged

ARKode with L:d:(Block)Vector #11660

Merged

sebproell mentioned this pull request Mar 17, 2021

Update SUNDIALS KINSOL interfaces to their current master #11908

Merged

bangerth mentioned this pull request Mar 31, 2021

Make KINSOL interfaces work for newer SUNDIALS versions #11829

Closed

bangerth mentioned this pull request Apr 20, 2021

Use N_Vector wrappers in KINSOL interfaces #12066

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Introduce SUNDIALS N_Vector module #11395

Introduce SUNDIALS N_Vector module #11395

sebproell commented Dec 18, 2020 •

edited

sebproell left a comment

sebproell Dec 18, 2020

sebproell Dec 18, 2020

sebproell Dec 18, 2020

sebproell Dec 18, 2020

sebproell Dec 18, 2020

sebproell commented Dec 21, 2020

peterrum left a comment •

edited

peterrum Jan 10, 2021

sebproell Jan 12, 2021

sebproell left a comment

sebproell Jan 12, 2021

sebproell Jan 12, 2021

sebproell Jan 12, 2021

sebproell Jan 12, 2021

sebproell Jan 12, 2021

sebproell Jan 21, 2021

sebproell Jan 12, 2021

sebproell commented Jan 12, 2021

sebproell commented Jan 21, 2021

luca-heltai commented Jan 22, 2021

luca-heltai commented Jan 22, 2021

sebproell commented Jan 22, 2021

		NVectorView<VectorType> n_vector(vector);
		auto id = N_VGetVectorID(n_vector);

		* @note N_VDestroy() should not be called on the resulting N_Vector since
		* this would lead to a double delete. Let the destructor do this work.

		auto vector = create_test_vector<VectorType>();
		const auto const_vector = create_test_vector<VectorType>();

		auto n_vector = make_nvector_view(vector);
		auto const_n_vector = make_nvector_view(const_vector);

Introduce SUNDIALS N_Vector module #11395

Introduce SUNDIALS N_Vector module #11395

Conversation

sebproell commented Dec 18, 2020 • edited

sebproell left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

sebproell commented Dec 21, 2020

peterrum left a comment • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

sebproell left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

sebproell commented Jan 12, 2021

sebproell commented Jan 21, 2021

luca-heltai commented Jan 22, 2021

luca-heltai commented Jan 22, 2021

sebproell commented Jan 22, 2021

sebproell commented Dec 18, 2020 •

edited

peterrum left a comment •

edited