[WIP] use SUNDIALS to compute reaction

[jdannberg]

First version of using SUNDIALS in the operator splitting routine. If you have any feedback on the general structure, let me know!
For now, I have ignored the possibility that temperature and composition might use different finite elements. I have also simply replaced the old version instead of checking if deal.ii is compiled with SUNDIALS.

I am also looking for feedback on how to pick the input parameters for ARKode.

• I have followed the instructions for indenting my code.

For new features/models or changes of existing features:

• I have tested my new feature locally to ensure it is correct.
• I have created a testcase for the new feature/benchmark in the tests/ directory.
• I have added a changelog entry in the doc/modules/changes directory that will inform other users of my change.

[bangerth]

Very cool. I think you're real close.

[bangerth]

I'd just drop this change -- this will minimize the difference in the output files.

[bangerth]

This is hard to read. It corresponds to this function:

      AdditionalData(
        // Initial parameters
        const double initial_time      = 0.0,
        const double final_time        = 1.0,
        const double initial_step_size = 1e-2,
        const double output_period     = 1e-1,
        // Running parameters
        const double       minimum_step_size                     = 1e-6,
        const unsigned int maximum_order                         = 5,
        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,
        // Error parameters
        const double absolute_tolerance = 1e-6,
        const double relative_tolerance = 1e-5);

Perhaps you can write this as

Suggested change

	SUNDIALS::ARKode<VectorType>::AdditionalData data (0.0, time_step, 0.001*time_step, 0.1*time_step, 1.e-6*time_step);
	SUNDIALS::ARKode<VectorType>::AdditionalData data (/* initial_time = */0.0,
	/* end_time = */ time_step,
	etc etc etc);

?

[bangerth]

Separately, though, it would avoid potential bugs if you replaced the initial time zero with the actual current time, and the end time as current time + time step.

I would probably choose a larger initial time step than 0.001time_step. Do you actually need any output every 0.1time step? This introduces a time point that ARKode needs to hit, and if you don't need the actual output, just set that parameter to zero or equal to time_step.

[jdannberg]

For the formatting: I had to place the comment at the end for the indentation to look reasonable.

[bangerth]

Move the declaration of this variable out of the loop over all cells. This way it needs to allocate memory only once.

[bangerth]

Same for this variable.

[bangerth]

It might be nice to create a lambda function of the form

  auto copy_fields_into_one_vector = [n_q_points,n_fields](const auto &in, auto &fields) {...};

that you place somewhere up above, and then you here just have to do

  copy_fields_into_one_vector (in_C, fields);

[bangerth]

This, too, could then be places into a local lambda function where then you just have to call copy_one_vector_into_fields(fields, in_C) or some such.

Names up for discussion of course.

[bangerth]

I don't recall if there was a way to ask the ARKode object for the number of steps it took. But if there isn't, then this will work:

unsigned int n_ode_timesteps = 0;
ode.output_step ([&n_ode_timesteps](const double       t,
                       const VectorType & sol,
                       const unsigned int step_number) {
  n_ode_timesteps = step_number;
});

[jdannberg]

I checked, and the solve_ode function returns the final number of computed steps. So I could use that, but it might be a different number for each cell. So the question is: What would be a reasonable number to return? The average? The maximum? The sum?

[jdannberg]

I addressed your comments (other than the last one about the steps), so my next question would be: how do we want to move forward with this? At the moment, we don't even have a section for the operator splitting in the timer output, so I do not know how this change affects the runtime. So that could be something to add in a separate PR?

We have some test cases with an analytical solution, so we do know that the new solution is more accurate now. But do we need to quantify this? Do we want to make the tolerance an input parameter?

Other than that, if this looks reasonable, I could now also add the implementation for the case where temperature and composition have different finite elements.

And finally, do we want to wait with this until we require SUNDIALS, or should I keep the old code and add a check if we have SUNDIALS?