Added two embedded methods

[brownbaerchen]

I added the first two embedded methods to give you a chance to complain before copy-pasting all the others.

The implementation has b.shape = (2, len(weights)). The weights property will give the first entry, which is the higher order method and the new weightsSecondary (notice the camel case...) will give b[1].
This means you can use embedded methods just like any other, totally ignoring the secondary method.

I test the secondary method by constructing a new RK scheme in the test with the b and order from the secondary method and the rest the same as the primary one.

Any more wishes before I add the rest of the methods?

[codecov-commenter]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 100.00%. Comparing base (4d4ab43) to head (4301f55).

Additional details and impacted files

@@            Coverage Diff            @@
##              main        #3   +/-   ##
=========================================
  Coverage   100.00%   100.00%           
=========================================
  Files           14        14           
  Lines          899       955   +56     
=========================================
+ Hits           899       955   +56     

Flag	Coverage Δ
smart-tests	100.00% <100.00%> (ø)

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[tlunet]

This looks quite fine on the technical side ... however looking at your changes I do wonder from the conceptual side :

from your implementation of Heun_Euler embedded (and others), it seems that you have the additional weights providing one order less (thought it was one order more, guess I had it wrong). So it means that any other
Q-generator could have those embedded weights, not only RK ? ... For instance : a collocation method could have naturally some embedded weights, that uses one accuracy order less than the quadrature formula. For RADAU-RIGHT and LOBATTO, then those embedded weights would simply be : [0, ...., 0, 1] !

If my interpretation is true, then it would be more generic to add this on the QGenerator level, rather than the RK one. Quick thinking on how to do that :

1. Add weightsEmbedded (or weightsSecondary, or else ...) property to the base class

class QGenerator(object):
    // ...
    @property
    def weightsEmbedded(self):
        raise NotImplementedError(f"no embedded weights implemented for {self.__class__.__name__}")

but don't impose its overriding in the register decorator.

2. In QGenerator.genCoeffs, add another argument like embedded=False that, when True, uses the weightsEmbedded property to extend the weights by the embedded weights, making it a 2darray of shape (2, M)

if embedded:
    out[1] = np.vstack(out[1], self.weightsEmbedded)

=> no need to check if embedded weights are implemented, we let the default parent method raise the NotImplementedError if not overriden

3. Add the embedded=True argument in the genQCoeffs function and pass it to the genCoeffs method

4. In RK, add a default bEmbedded (or b2, or else) class attribute that is default to None, and transformed into a numpy array in the checkAndStore decorator only if it not None. Then override the weightsEmbedded property that raises NotImplementedError if bEmbedded is None, else return those.

Then, the embedded order would simply be obtained as QGenerator.order - 1, which may not require a dedicated property after all ...

And for testing (in test_base or test_convergence), run the tests on all $Q$-generators, but wrap those into a try catch block expecting NotImplementedError, which will basically pass. In particular, consistency tests (correct size, etc) can simply be added in testAdditionalCoeffs, while convergence tests simply need to be added in testDahlquist, where an additional useEmbedded=False argument should be added to the errorDahlquist method (and probably to the solveDahlquist too), all of this wrapped in a try catch block expecting NotImplementedError.

Of course : works only if my first interpretation is true 😅 ... But if yes, then this would not modify a lot your approach, just make it way more generic and avoid differentiating in butcher.py between embedded and non-embedded classes. What do you think ?

[brownbaerchen]

That's a good idea, actually!

[tlunet]

Then there is maybe also a way to add your work on embedded SDC methods in an advanced tutorial 😉 ...

[tlunet]

Of course, idea of later 😅 !

[tlunet]

Once this is finalized and all RK schemes from pySDC are added, then I'll create the first beta release 0.1.0, cf #6

[brownbaerchen]

We need the property orderSecondary because the difference in order can be larger than one as is the case in the additional method I added.

I followed the steps you suggested, but for some reason the order is only 1 with b = [0, ..., 0, 1]. Don't know why and I don't really have bandwidth for figuring this out right now. Maybe you have an idea how to fix this?

[tlunet]

There is probably not need to implement the embedded Q yet : that was just a hunch on my side, and probably need to be more analyzed before implementing it ... I think it make sense for SDC, but probably not for direct collocation, since the collocation with or without the quadrature update provide the same solution for Lobatto and Radau-Right.

So you can probably leave the default non-implemented embedding for Collocation, and juste use if for RK methods for now.

PS : there is probably more to it, actually ... see embedded coefficients for Lobatto I to III. In fact, I would even see it as a student project proposal, taking some of the content from this project (generic discontinuous collocation methods), and combining this with generic embedding (and SDC eventually ...)

[tlunet]

Finally, about the orderSecondary (or orderEmbedded ...) : I think we could implement a default property in QGenerator that return self.order-1, and only override when embedded order is not order-1 (seems like it corresponds to special cases ...).

[tlunet]

you may wanna merge main into your branch to avoid the tutorial test error ... the PR will be squashed anyway

[tlunet]

Hum, that's the weirdest error I've ever seen ... nvm, I'll merge to main anyway and we'll solve that after.

What about secondary VS embedded ? (sry, forgot that I had to apply the review to make it visible to you 😅 )

[tlunet]

Looks fine enough, I'll correct in additional commit later ...

Thanks a lot for the work !