Refactor of update_nodes

[brownbaerchen]

I started reading a book called Clean Code and decided to give some of the things I learned a shot. I refactored the update_nodes functions in the generic_implicit and imex_1st_order sweepers. Feel free to reject this. What constitutes good code is subjective after all and I can see how someone could find the current version more readable than what I suggest.

The principles I tried to apply are mainly:

• Forgo comments in favour of expressive function and variable names
• Make small functions that do as little as possible

A benefit gained from the abstraction of small functions is that the imex sweeper can be shortened significantly. Of course, this comes at the expense of difficulty to understand abstract code. You cannot really understand what the imex sweeper does without considering the above level of abstraction. On the other hand, the differences to the generic implicit sweeper become much more clear as you don't really need to implement any shared functionality.
If you like these changes, I would apply them to some of the other sweepers as well and I believe that the explicit sweeper or the RK sweeper can be shortened significantly as well. In this case, it would be worth considering moving some changes to the core sweeper module.

At first I didn't want to start refactoring existing code. However, I may have to start taking memory a bit more seriously and this is a side effect of these changes. Memory usage is indeed reduced because the current implementation generates a right hand side for the problem to solve from a copy of the "integral" at the respective node. In this PR, I want to change that to construct a right hand side vector from the integral and so on and directly pass this to the solver.
Another side effect is gained efficiency because I suggest to compute $(Q-Q_\Delta)F$ instead of $QF - Q_\Delta F$ as in the current implementation. I actually have a separate copy of these sweepers in my project folder right now which implement this already. By refactoring some other stuff a bit more, I could get rid of them and we could have the more efficient implementation as default.

Please do comment! Do you like this style of coding or do you like a different level of abstraction that generates more self-contained files? Do you require comments and docstrings? Again, I don't claim to have objectively improved the prettiness of the code in this PR.

While we're on the topic of more expressive naming, I suggest to rename a few things, for instance

• integrate -> integrate_right_hand_side
• solve_system -> Euler_step
• description -> configuration
• remove "classes" in the path to implementations, i.e. pySDC.implementations.sweepers.generic_implicit etc.
• ...
  How do you feel about such things?

[tlunet]

I like the initiative 😄

I would avoid the snake_case tho (according to the current coding guidelines ... 😉).

Agree on removing the "_classes", would also remove underscores and upper cases in any module classes, simplify the names for the sweeper module (e.g implicit instead of generic_implicit), put the name of the classes in PascalCase, and import all implemented classes in the pySDC.implementations.sweepers.__init__.py, such that one could simply use them like this (for instance) :

from pySDC.implementations.sweepers import Implicit, IMEX, MultiImplicit

Also, I would not specify the ODE order for the sweeper when first order, since it is how SDC is implemented by default ...

[pancetta]

I agree with @tlunet. This would remove quite a few poor initial-stage decisions, indeed. Further comments:

• Unless you deviate too much from the original formulas found in the papers, feel free to change the sweepers to make them more efficient.
• This is still a big playground, not a high-performance codes.. people should be able to understand the code "easily" (maybe with the controllers as an exception). Don't abstract everything away, if people do not recognize the actual algorithms anymore.
• Don't call it Euler_step this is confusing (why Euler).
• One thing we could also consider is that the residual, once computed, can be used to convergence AND for the sweeper. Currently we don't exploit that and technically this may be difficult to do well, but it's an option.

[brownbaerchen]

I agree with @tlunet. This would remove quite a few poor initial-stage decisions, indeed. Further comments:

• Unless you deviate too much from the original formulas found in the papers, feel free to change the sweepers to make them more efficient.
• This is still a big playground, not a high-performance codes.. people should be able to understand the code "easily" (maybe with the controllers as an exception). Don't abstract everything away, if people do not recognize the actual algorithms anymore.
• Don't call it Euler_step this is confusing (why Euler).
• One thing we could also consider is that the residual, once computed, can be used to convergence AND for the sweeper. Currently we don't exploit that and technically this may be difficult to do well, but it's an option.

Actually, the abstraction is meant to aid seeing the algorithm in the code. So, for instance, instead of

# gather all terms which are known already (e.g. from the previous iteration)
# this corresponds to u0 + QF(u^k) - QdF(u^k) + tau

# get QF(u^k)
integral = self.integrate()
for m in range(M):
    # get -QdF(u^k)_m
    for j in range(1, M + 1):
        integral[m] -= L.dt * self.QI[m + 1, j] * L.f[j]
        # add initial value
        integral[m] += L.u[0]
        # add tau if associated
        if L.tau[m] is not None:
            integral[m] += L.tau[m]

I want to call a function self.build_right_hand_side(), which has the following code

def build_right_hand_side(self):
     rhs = self.initialize_right_hand_side_buffer()
     self.add_Q_minus_QD_times_F(rhs)
     self.add_initial_conditions(rhs)
     self.add_tau_correction(rhs)
     return rhs

The idea is that you can actually read the algorithm in the code! Perhaps an even better example is the add_new_information_from_forward_substitution function. The downside being that the actual implementation of some of these functions is quite far removed in the file from the more high level things, requiring many jumps. On the other hand, we can put a lot of stuff in the core sweeper. Afaik, the update_nodes function always has these steps:

• build right hand side from last iteration
  • add $(Q-Q_\Delta)f$
  • add $u_0$
  • add tau correction
• sweep
  • forward substitution
  • solve
  • evaluate $f$

The differences, on the other hand are only in adding $(Q-Q_\Delta)f$, forward substitution and the solve step. So maybe we can write more shared structure into the core sweeper module and only overload the things that are actually different.
There are some exceptions of course. For instance, diagonal SDC sweepers should avoid calling the forward substitution function for performance reasons. But I think it is sensible to generate default behaviour and overload the default behaviour in edge cases rather than define no default behaviour and copy-paste a bunch of code.

Concerning renaming solve_system to Euler_step, my reasoning was that, to the best of my knowledge, we always solve implicit Euler, explicit Euler or IMEX Euler steps when calling this function. I guess it is possible to replace this with some other RK method, so maybe we should leave it more general, you're right.

[brownbaerchen]

Regarding computing the residual only once, I actually vote against it. It assume it is more efficient to compute $(Q-Q_\Delta)F$ and then to compute the residual only on the last node. Or not at all, if the algorithm doesn't require it. If we want to reuse $QF$, we have to separately subtract $Q_\Delta F$, I don't think we can gain much.

[lisawim]

• I agree with all of you! I also would prefer using the PascalCase due to the coding guidelines @tlunet mentioned (I already had this longer in my mind that we should change that. Now, it could be a good possibility to do that.).
• The current structure in the sweepers was very helpful for me to understand what's going on there. For better readability it seems to be good practice to move doings in functions. Thus I really like these ones such as add_Q_minus_QD_times_F(rhs), self.add_initial_conditions(rhs), self.add_tau_correction(rhs) and so on (although these are only a few lines, but several lines can be saved here when implementing new sweepers). These function names also indicate what's going on there, and I like that as well! From the SDC/PFASST papers it should be already intuitively clear what these functions do. But only for folks really dealing with that stuff. For all others there is one solution to make sure that the user understands how sweepers work in pySDC: writing a very detailed documentation! However, this documentation has to be written only once for the core sweeper class.
• Regarding the renaming solve_system -> Euler_step: I'm of two minds about that. You have your right-hand side consisting of the initial condition the parts where the integration is done, and the tau correction. And then you solve an (implicit) system to get the value for the next node. This could be indeed confusing for some. Otherwise, we can think about an implicit Euler. What we are doing here? Consider the implicit Euler of the form $$u_{n+1} = u_n + \Delta t f(u_{n+1})$$ for some problem $u'=f(u)$. In order to find the value for the next time step, the (probably nonlinear) problem $$u_{n+1} - \Delta t f(u_{n+1})=u_n.$$ has to be solved, which can be done by Newton, for instance. Now, when we think of what's going on in SDC, we also want to solve a system like $$u_{m+1}^{k+1} - \Delta t q_{m+1,m+1}^{\Delta}f(u_{m+1}^{k+1}) = u_n + \Delta t \sum_{j=1}^m q_{m+1,j}^{\Delta}f(u_j^{k+1}) + \Delta t \sum_{j=1}^M (q_{m+1,j} - q^{\Delta}_{m+1,j})f(u_j^k),$$ where $u_n$ is the initial condition at time $t_n$.

[lisawim]

So, this could be a similarity? In general, a system $$u_{n+1} - \Delta t c f(u_{n+1})=b$$ is solved. For implicit Euler, the system to be solved can be defined by $c=1$ and $b=u_n$. For SDC, the system so be solved is defined by $c=q_{m+1,m+1}^\Delta$ and $$b=u_n + \Delta t \sum_{j=1}^m q_{m+1,j}^{\Delta}f(u_j^{k+1}) + \Delta t \sum_{j=1}^M (q_{m+1,j} - q^{\Delta}_{m+1,j})f(u_j^k).$$ Does that make sense in some way? Anyway, this is what I would see when I read Euler_step..

[brownbaerchen]

So, this could be a similarity? In general, a system un+1−Δtcf(un+1)=b is solved. For implicit Euler, the system to be solved can be defined by c=1 and b=un. For SDC, the system so be solved is defined by c=qm+1,m+1Δ and b=un+Δt∑j=1mqm+1,jΔf(ujk+1)+Δt∑j=1M(qm+1,j−qm+1,jΔ)f(ujk). Does that make sense in some way? Anyway, this is what I would see when I read Euler_step..

This was exactly my reasoning for suggesting the rename. I understand now that Euler is perhaps not general enough, however solve_system is a bit too unspecific for my liking. But we don't need to fix that now. If someone comes up with a good idea, we can change it. Otherwise, we can just leave it as is for now.

[pancetta]

So, this could be a similarity? In general, a system un+1−Δtcf(un+1)=b is solved. For implicit Euler, the system to be solved can be defined by c=1 and b=un. For SDC, the system so be solved is defined by c=qm+1,m+1Δ and b=un+Δt∑j=1mqm+1,jΔf(ujk+1)+Δt∑j=1M(qm+1,j−qm+1,jΔ)f(ujk). Does that make sense in some way? Anyway, this is what I would see when I read Euler_step..

This was exactly my reasoning for suggesting the rename. I understand now that Euler is perhaps not general enough, however solve_system is a bit too unspecific for my liking. But we don't need to fix that now. If someone comes up with a good idea, we can change it. Otherwise, we can just leave it as is for now.

Well, the only purpose of this function is to "solve a system". This does not have a focus on time-stepping (at least from the perspective of the API), but on linear or nonlinear solvers. This is why it is called solve_system.

[brownbaerchen]

Ok, let's forget about renaming for a moment and focus on this PR. :D @pancetta, how do you feel about it? How do you feel about moving some of the changes to the core module? I am not offended if you don't feel favourably :D

[pancetta]

I like the renaming of things and a certain level of abstraction, too, but I'm not sure how this would look like camelCase-style, esp. those long descriptive function names. What do you think @tlunet, @brownbaerchen?

[brownbaerchen]

Do you also prefer underscores? I thought I was the only one at this point :D But I will give in to the camel case crowd if needed. I can see some benefit in that...

[tlunet]

Not really sure about this thing :

def build_right_hand_side(self):
     rhs = self.initialize_right_hand_side_buffer()
     self.add_Q_minus_QD_times_F(rhs)
     self.add_initial_conditions(rhs)
     self.add_tau_correction(rhs)
     return rhs

When each method uses basically one or two line of code, it's probably better to have a comment in the code and then the line. I think the previous version is most understandable

[pancetta]

Do you also prefer underscores? I thought I was the only one at this point :D But I will give in to the camel case crowd if needed. I can see some benefit in that...

I guess I do prefer underscores (as you can see all over pySDC), but I understand @tlunet's point in the contributors guide. I don't care that much, but for long function names underscores might be better to read.

[tlunet]

Concerning the solve_system or euler_step, it always solves something like this :

$$ u - \alpha F(u, t) = rhs $$

which is the base solver for any Backward Euler step, or any implicit SDC update, etc ... Since it's linked to the problem itself (and it's f function), why not simply something like this :

    def solve(self, alpha, b, t):
        r"""
        Solve the following (non-) linear system :

        .. math::
            u - \alpha F(u, t) = b

        with :math:`F(u,t)` the RHS function of the problem.

        Parameters
        ----------
        alpha : float
            The :math:`\alpha` coefficient
        b : dtype_u
            Right-hand side :math:`b` for the system
        t : float
            Time of the solution.

        Returns
        -------
        u : dtype_u
            Solution of the (non-)linear system
        """
        raise NotImplementedError() 

Then it can be used for any implicit updated in pySDC. and also as for the generic exact solution implemented by @brownbaerchen, we could have a generic implementation of this method in the base Problem class using some generic non-linear solver of scipy (that could switch to a newton_krylov based solver as soon as the number of degree of freedoms gets higher than a given threshold.

PS : on a side note, I agree that long function name are not great in camelCase, but I also believe that long function name are not great either. So forcing camelCase also forces to find names that are not too long, not too short (which is not very enforced when we allow ourselves snake_case ...)

Also, nothing prevents to add a short line of documentation in the header of the function, ex :

def doSomething(self):
    """Do the something we talked about"""

[pancetta]

Yes, I agree! Both with the solve and the naming scheme. I know, Clean Code suggests otherwise, but hey, they are doing SE, not RSE.

[brownbaerchen]

Since my spider sense picks up only limited enthusiasm about this style of coding, I will stop with the refactoring for now ;)