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

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.

-
@gnzlbg: Perhaps you could include a reference or two to the methods you're discussing? This is the sort of question our community tends to be good at, so the lack of responses indicates that we're not sure exactly what you're asking. – Aron Ahmadia Aug 3 '12 at 12:01
@AronAhmadia done. – gnzlbg Aug 5 '12 at 13:07
Thanks for the added information. I don't have an answer, but I believe that the objections to the method of lines on the page you've linked to are misguided. They are actually just objections to the use of multistep time discretizations in combination with dynamically adaptive spatial grids. Those choices are independent of whether one uses the method of lines or the method of Rothe. Within the method of lines, one can obtain high order, one-step schemes by using Runge-Kutta time discretization. – David Ketcheson Aug 5 '12 at 13:49
@DavidKetcheson Yes, I dont agree completely with the arguments given in the discussion either. Still, the method of Rothe really is better for FEM discretizations of convection dominated problems: the stabilization parameter depends on both the time step, the spatial, and the temporal discretization used. Using the method of lines your stabilization parameter would get tricky to say the least. However, pros and cons of each method is a whole topic in itself, and one would have to consider space-time formulations too. – gnzlbg Aug 5 '12 at 14:03
Since @Wolfgang wrote the Method of Rothe description you linked to as part of documentation for deal.II, I've asked him to try and tackle your answer. – Aron Ahmadia Aug 7 '12 at 15:06
show 1 more comment

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.

-
This "answer" is more of a comment -- I don't see how it answers the question. Also, in the method of lines, if you use a one-step method, there's nothing wrong with having a different set of ODEs at every step. – David Ketcheson Aug 8 '12 at 6:45
@WolfgangBangerth "One-step" method as opposed to multi-step method. Runge-Kutta methods, including Rosenbrocks, DG-in-time, and IMEX variants are easy to use with adaptivity. Design of time integration schemes can be similarly subtle to spatial discretization, considering strong stability properties, embedded error control where the embedded schemes have desired stability properties, adjoints with desired properties, targeting specific parts of the spectrum, obtaining a desired stage order, preserving geometric properties, being as inexpensive for the implicit solver as possible, etc. – Jed Brown Aug 8 '12 at 7:37
You can of course do this. The problem with your approach is that you need to interpolate (or project somehow) solutions on previous meshes onto the the new mesh. This introduces an error that is not just a Galerkin projection and is, thus, not easy to analyze -- the result is not simply a Galerkin scheme. – Wolfgang Bangerth Aug 8 '12 at 21:18
On the other hand, if you derive a method via the Rothe approach, you get a PDE for each time step in which the previous time step's solution appears on the right hand side. By discretizing, you test this previous solution with test functions on the current mesh. This can be interpreted as some sort of projection -- but it is not the L2 projection but instead something that is a weighted combination of L2 and H1 projection with the weight depending on the time step. Since it's a pure Galerkin scheme, it is easy to analyze and I would venture the guess that the error is smaller this way. – Wolfgang Bangerth Aug 8 '12 at 21:20
Now, as a final point, my best guess is that if you only change the mesh occasionally, it's not going to make much of a difference in practice. Since testing a solution on one mesh with test functions from a different mesh is sort of a pain, I do this interpolation all the time -- without too many problems ;-) – Wolfgang Bangerth Aug 8 '12 at 21:21