Jacobian-free Newton-Krylov implicit solver for diagonally-implicit RK schemes

[cpethrick]

PR for introducing implicit RK using Jacobian-free Newton-Krylov (JFNK) for the implicit solver. JFNK can be quite fast as no Jacobian matrix needs to be stored. Rather, the Jacobian-vector product which arises during Newton iterations is approximated on-the-fly using finite differences. The absence of a Jacobian implies that automatic differentiation is not required. In the scope of this PR, there will be no preconditioner used.

This PR has also improved the RK ODE solver to use better object-oriented programming practices.

Changes related to JFNK:

• Rename ExplicitODESolver to RungeKuttaODESolver, as both implicit and explicit RK are handled by the same function.
• Add diagonally implicit RK (DIRK) to RungeKuttaODESolver
• Create a JFNK solver using deal.ii's GMRES solver. Currently this and all associated classes are in src/ode_solver/JFNK_solver. JFNK_solver.<cpp/h> holds deal.ii's GMRES solver. jacobian_vector_product.<cpp/h> exists mainly to act as a matrix-like operator for the GMRES solver, and also computes the unsteady residual.
• Added tests to demonstrate implicit RK convergence order for 1D advection and burgers'.

Other changes to clean up RungeKuttaODESolver:

• Added parameters for JFNK to LinearSolverParameters, re-using GMRES parameters where possible.
• RungeKuttaODESolver and derived class RRKExplicitODESolver now assign n_rk_stages as a template argument for better performance. create_RungeKuttaODESolver has been added to ODESolverFactory has been modified to assign n_rk_stages based on parameters.
• Changed from selecting an RK method by its order to named methods using the parameter runge_kutta_method. This allows multiple RK methods to have the same order, and streamlines the addition of methods.
• The assignment of the Butcher tableau has been moved to classes derived from RKTableauBase. An associated factory, create_RKTableau, has been made.

[cpethrick]

I added the "c" part of the butcher tableau to all of the RK methods. This is unused in the current RK solver, but I have added it because PR #111 uses it for some unsteady manufactured solutions.

[cpethrick]

In the development of this PR, I also tested implicit RK using the AD-based linear solver already implemented in PHiLiP. I thought it would be interesting to compare the two now that JFNK is functioning well. I modified runge_kutta_ode_solver.cpp (see attachment) such that I could compare AD to JFNK for the test 1D_TIME_REFINEMENT_STUDY_ADVECTION_IMPLICIT. A summary of the comparison is as follows:

• I needed to decrease the size of the first timestep from 0.2 to 0.01 in order for AD to converge. So, for the test set-up, JFNK is more robust. As a caveat, the timestep for AD to converge could likely be increased by increasing the GMRES tolerance.
• JFNK and AD both converge at the theoretical convergence rate, shown in the attachments.
• JFNK has an L2 error about 6 times smaller than AD.
• JFNK takes more than 2x the time taken by AD. This is to be expected for a serial computation - JFNK uses many more operations than AD (which can be faster on GPU calculations). However, the efficiency hasn't been optimized for either method.

I've attached the modified runge_kutta_ode_solver.cpp which uses AD, the .prm file with reduced initial timestep, as well as the full output of 1D_TIME_REFINEMENT_STUDY_ADVECTION_IMPLICIT using AD and JFNK.

runge_kutta_ode_solver_ADversion.cpp.txt time_refinement_study_advection_implicit.prm.txt time_refinement_study_advection_implicit_AD_convergence.txt time_refinement_study_advection_implicit_JFNK_convergence.txt

[dougshidong]

Nice analysis.

As a side note, you can assemble the matrix-vector product using AD such that it becomes a JFNM. The real comparison is AD vs Fréchet FD.

As you mentioned both have opportunities for code optimization to speed them up. The real saving we're obtaining is the memory usage as we get to higher and higher orders.

[cpethrick]

I was discussing with @jbrillon how to best set up this PR for compatibility with PR #111. We decided that it would be easiest to rename the files runge_kutta_ode_solver.<cpp/h> back to explicit_ode_solver.<cpp/h> such that GIT recognizes it as the same file. The parameter names, etc. are unchanged from previous commits. I have opened an issue (#174) to change the names to runge_kutta_ode_solver.<cpp/h> once #111 has been merged into master.

[dougshidong]

What's left for this PR to be merged?

[cpethrick]

@dougshidong It should be ready to merge soon (within a week, I would guess). @jbrillon was going to do a final review in the coming days. I have confirmed that ctest behaves the same as before, and the time taken for explicit timestepping has not changed.