PDE discretization with the method of rothe and the method of lines (Modular implementation)

Question

The Heat equation is discretized in space with FV (or FEM), and a semi-discrete equation is obtained (system of ODEs). This approach, known as the method of lines, allows to easily switch from one temporal discretization to another, without code duplication. In particular, you can reuse any time integrator for ODEs without much effort. This is very convenient because if you decide to change your spatial discretization from FV to say FE, you still get a semi-discrete equation and your time integrators still work.

Now I am trying to implement the method of rothe for the same problem. However, discretizing in time first forces me to rewrite the spatial discretization for every temporal discretization scheme I might want to use. This eliminates the reuse of time integrators that I previously had, and makes it very complicated to write modular software that can discretize a PDE using both the method of lines or the method of Rothe.

Is there a way of implementing both approaches, without code duplication?

Edit:

In convection dominated problems, the FE discretization needs stabilization both in time and space, making the method of Rothe the "best" choice. However, this is not the case for FV/DG methods.

In the method of lines, the PDE is discretized first in space, and then in time. In the method of Rothe, the PDE is discretized first in time, and then in space. The third possibility is to discretize both in space and time simultaneously (also known as space-time discretizations). A discussion about the method of lines and the method of Rothe can be found here. For more information the book "Finite Element Methods for Flow Problems" from Donea and Huerta is a good resource.

Answer 1

I don't really have much more to say than I already did on the pages you linked to, but for me the primary arguments go like this:

• In many problems, one needs to adapt the mesh between time steps. The conceptual framework for doing this is the Rothe method where one can choose the spatial discretization independently at each time step whereas the method of lines a priori assumes that the PDE is converted into a system of ODEs -- which is incompatible with mesh adaptation.

• On the other hand, if you don't want to adapt your mesh, then it doesn't matter: in most cases, if the spatial discretization is going to be the same between time steps, then it doesn't matter whether you want to first discretize in space and then in time using your favorite time integrator, or the other way around -- you will come out with the same discrete problem to solve in every time step. In such cases, the Rothe method and the method of lines are then the same.

• This extends to the case where you only want to adapt the mesh every once in a while: you can consider this as the method of lines applied to a number of time steps, then adapting the mesh, then one more set of time steps where you apply the method of lines. Or you can think of this as the Rothe method where you just happen to only adapt the mesh every once in a while. It will essentially come out to the same numerical scheme, just a different philosophical viewpoint.

It might be worthwhile adding one more point: In the ODE world, one often uses high order schemes with several stages or multiple steps. Thus, there is significant benefit to bottling up these algorithms into packages that you only have to hand an ODE system in one way or another. On the other hand, for time dependent PDEs, most of the time one uses rather simple time stepping methods (with the notable exception of some hyperbolic solvers): for example, Crank-Nicolson, BDF-2 or just the forward or backward Euler schemes. For these simple time integrators, it isn't particularly difficult to hand-code the time integration since it is so much simpler than the spatial discretization. What this means is that the price to pay for thinking in terms of the Rothe method -- not being able to use an ODE solver package -- is a small one, whereas the price to pay for using the method of lines -- not being able to adapt the mesh between time steps -- is a large one. That might explain why most people in the adaptive finite element world prefer to think along the Rothe method.

As a corollary to the last point, and coming back to the original question: It is indeed true that it is difficult in the Rothe method to package everything up nicely in an object oriented way. However, (i) as long as you stick with a single class of ODE integrators, you can of course still tabulate the coefficients of the various stages of the integrator in a class and have the code that computes them be entirely encapsulated; and (ii) the fact that one commonly uses relatively simple time integrators for time dependent PDEs means that the amount of code necessary to implement the time integration is usually vastly smaller than the amount of code that deals with the spatial discretization. In other words, I know of no way how to nicely separate the spatial from the temporal discretization using the Rothe method, but I haven't usually found this to be too much of an obstacle since so little code is necessary to write, say, a BDF-2 or RK-4 scheme.