Creating and destroying work arrays in SEULEX solver

[dholland88]

When using the SEULEX ODE solver in the cfd plugin, there is also a significant amount of work being done by creating and destroying the work vectors every time the equations are being solved. It would be give a decent performance improvement if this could be refactored so that they are only created once and destroyed once at the end of the CFD simulation. This would just require that the same species are being solved for in every cell for every time step, which is generally true.

I have not checked/tested the CVODE solver for a similar issue.

[whitesides1]

Thanks for the report. Do you have an estimate for the amount of time spent in these allocs/deallocs relative to the solve (including mechanism size and platform)?

In general, you are right on the requirement, but the refactoring effort would be significant, based on how the solver classes are designed.

[dholland88]

From the profiling I did it spends about ~10% of its time doing allocs/deallocs.

Yes, I did look at doing it myself but soon realised it was going to be a bigger job than I had hoped.

[whitesides1]

OK, I've thought about this a bit and don't have a nice way of refactoring this and keeping the flexibility of the current structure (especially when considering GPU solvers).

In my test cases, it didn't seem like it was it was as much as 10%, so if you can provide a reproducer for the conditions where you see the highest impact of memory allocation, I can use it if/when there is a larger push to refactor.

If you're interested you can try the following patch. In my testing it didn't make a big difference, but maybe it will in your case.

diff --git a/applications/cfd_plugin/interfaces/seulex_cpp/seulex_cpp.cpp b/applications/cfd_plugin/interfaces/seulex_cpp/seulex_cpp.cpp
index 993c901..321668e 100644
--- a/applications/cfd_plugin/interfaces/seulex_cpp/seulex_cpp.cpp
+++ b/applications/cfd_plugin/interfaces/seulex_cpp/seulex_cpp.cpp
@@ -627,22 +627,51 @@ void seulex::setup_step_sequence(std::vector<int>& nj, std::vector<double>& a)
 
 void seulex::create_work_vectors(N_Vector y)
 {
-  yh = N_VClone(y);
-  dy = N_VClone(y);
-  dyh = N_VClone(y);
-  del = N_VClone(y);
-  wh = N_VClone(y);
-  scal = N_VClone(y);
-  tmp1 = N_VClone(y);
-  t.resize(km+1);
-  for(int i = 0; i <= km; ++i)
-  {
-    t[i] = N_VClone(y);
+  if(N_VGetVectorID(y) == SUNDIALS_NVEC_SERIAL) {
+    int num_elems = N_VGetLength(y);
+    double* vector_memory = (double*) malloc(num_elems*(7+km+1)*sizeof(double));
+    yh = N_VCloneEmpty(y);
+    N_VSetArrayPointer(vector_memory, yh);
+    dy = N_VCloneEmpty(y);
+    N_VSetArrayPointer(vector_memory+1*num_elems, dy);
+    dyh = N_VCloneEmpty(y);
+    N_VSetArrayPointer(vector_memory+2*num_elems, dyh);
+    del = N_VCloneEmpty(y);
+    N_VSetArrayPointer(vector_memory+3*num_elems, del);
+    wh = N_VCloneEmpty(y);
+    N_VSetArrayPointer(vector_memory+4*num_elems, wh);
+    scal = N_VCloneEmpty(y);
+    N_VSetArrayPointer(vector_memory+5*num_elems, scal);
+    tmp1 = N_VCloneEmpty(y);
+    N_VSetArrayPointer(vector_memory+6*num_elems, tmp1);
+    t.resize(km+1);
+    for(int i = 0; i <= km; ++i)
+    {
+      t[i] = N_VCloneEmpty(y);
+      N_VSetArrayPointer(vector_memory+(7+i)*num_elems, t[i]);
+    }
+  } else {
+    yh = N_VClone(y);
+    dy = N_VClone(y);
+    dyh = N_VClone(y);
+    del = N_VClone(y);
+    wh = N_VClone(y);
+    scal = N_VClone(y);
+    tmp1 = N_VClone(y);
+    t.resize(km+1);
+    for(int i = 0; i <= km; ++i)
+    {
+      t[i] = N_VClone(y);
+    }
   }
 }
 
 void seulex::destroy_work_vectors()
 {
+  if(N_VGetVectorID(yh) == SUNDIALS_NVEC_SERIAL) {
+    double* vector_memory = N_VGetArrayPointer(yh);
+    free(vector_memory);
+  }
   N_VDestroy(yh);
   N_VDestroy(dy);
   N_VDestroy(dyh);