[WIP] Quadrature Revamping #86

torbjoernk · 2014-09-04T07:35:08Z

Here is a proposed layout for the new quadrature. Only the quadrature tests are supposed to compile. The sweepers and examples have yet to be altered to not use "augmented nodes" any more.

The Q and S matrices are now square and nodes only holds "real" quadrature nodes. The quadrature_traits template is used to specify whether the left/right node are the start and/or end of the time interval. This is in preparation for #71 and #31.

To shorten the source files, I've split up the different quadrature types into separate files in a separate namespace.

splitting into separate files in subdirectory note: just include 'quadrature.hpp' for all the quadrature stuff

memmett · 2014-09-04T17:15:54Z

I like where this is going but need to be educated a bit...

Can you give an example of how the trait mechanism for left_is_node and right_is_node works, and contrast it to using the left_is_node and right_is_node methods? Given that the sweepers should be dynamic, is there any situation (internal to PFASST) where the trait mechanism would be used instead of the methods?

Splitting up all the different types of quadrature seems to have resulted in a boat load of boiler plate code -- can this be shortened up significantly?

Is there some fundamental reason why the matrices are square (aside from making it easier to compute Q*Q)?

Do we really need those m_ prefixes?

We need to fix that UNUSED(to), but that can wait...

pancetta · 2014-09-04T18:54:22Z

No fundamental reason for Q to be quadratic, but IMHO much more consistent and easier to handle (e.g. for Q*Q or St. Martin's trick). Plus: no special treatment (in Q and number of nodes) for quadrature rules where the left boundary point is not a quadrature node.

mlminion · 2014-09-04T19:11:51Z

How do you propose to handle Guass-Legendre nodes with a square matrix
the size of the number of nodes? You would need an extra quadrature rule
to integrate to t_{n+1}.

M

On 9/4/2014 11:54 AM, Robert Speck wrote:

No fundamental reason for Q to be quadratic, but IMHO much more
consistent and easier to handle (e.g. for Q*Q or St. Martin's trick).
Plus: no special treatment (in Q and number of nodes) for quadrature
rules where the left boundary point is not a quadrature node.

—
Reply to this email directly or view it on GitHub
#86 (comment).

pancetta · 2014-09-04T20:09:57Z

Quite true, but we agreed on a special treatment for the boundaries anyway. We have a vector q which includes the weights from t_n to t_{n+1} and is used only if the last node is not equal to t_{n+1}. Similarly, we test whether the first node is equal to t_n and if not, we use the first row of Q (otherwise, we don't, so no adding up zeros).

The reason for that: SDC is an iterator for u at the nodes tau_m only. u_1 = u(t_n) and/or u_M = u(t_{n+1}) are special cases, where we can save a few operations and we do so using left_is_node and right_is_node for convenience.

memmett · 2014-09-04T20:13:46Z

@pancetta Can you just clarify: if the user asks for 3 Gauss-Legendre nodes, do we get a 3x3 matrix?

pancetta · 2014-09-04T20:17:01Z

@memmett exactly, you will always get an MxM matrix Q, an MxM matrix S and the weights vector q for the full interval , no matter what kind of nodes you request. First row might be all zeros in S and Q, but who cares.. we can (and will) catch those cases in the sweeper.

And the order of your quadrature is always directly linked to M, i.e. 2M for Gauss-Legendre, 2M-2 for Gauss-Lobatto. Principle of least surprise :-).

memmett · 2014-09-04T20:26:07Z

Yep, LGTM except for those technical questions that I raised...

mlminion · 2014-09-04T20:26:37Z

Oh please don't start calling the left endpoint u_1!

u_1 = u(t_n)

pancetta · 2014-09-04T20:27:42Z

@mlminion u_1 was meant to be the first node... sorry for the confusion.

mlminion · 2014-09-04T20:35:22Z

On 09/04/2014 01:27 PM, Robert Speck wrote:

@mlminion https://github.com/mlminion u_1 was meant to be the first
node... sorry for the confusion.

—
Reply to this email directly or view it on GitHub
#86 (comment).

The first node is very often at t_n, so I am not sure what you mean.

I like the notation in the MLSDC paper (see the footnote at the bottom
of page 3 if you are confused). You always have M+1 nodes and M
substeps. t_0=t_n and t_M=t_{n+1}.

I don't understand how having matrices of different sizes is going to
simplify anything.
The sweepers can do the right thing when rows or columns are zero, and you
don't have to drag around an extra quadrature rule because it is always
there.

M

torbjoernk · 2014-09-05T06:57:30Z

Thanks for your remarks. I'll comment on the traits and boilerplate code.

The traits are meant only for internal use and were the result of my inability to achieve the same effect with plain (static) const bool variables left_is_node and right_is_node. (See this Gist)
The problem I had was the difference in accessing the static class variable VAL and the overwritten method value().
When I rethink this now, I guess, it's a little easier to read without the traits (as we need to override the two methods anyway). 😉 I'll give it a try.

Regarding the boilerplate code, you probably mean all the copy and move constructors and operator= overloading, right? I'll try to remove and avoid that. Specializations of these is probably not needed as the quadrature classes do not handle any resources, thus the default (compiler-generated) implementations should be sufficient.

using naming convention of rest of codebase

why does clang not complain? o.O

everything should still work as before

also adjusted AdVec Vanilla SDC example note: only vanilla_sdc compiles and runs

memmett · 2014-09-05T19:25:43Z

include/pfasst/interfaces.hpp

@@ -150,6 +154,12 @@ namespace pfasst
       */
      virtual void advance() = 0;

+      virtual void do_inner_nodes(const size_t m, const time t, const time ds) = 0;


Do these belong at the interface level?

memmett · 2014-09-05T19:37:13Z

@torbjoernk The changes you made to remove traits and most of the boiler plate code look good. I made a few comments inline.

@mlminion @pancetta I think there is a conservation-of-pain law lurking around here: with either convention (left/right always included vs not) the sweeper will have some extra logic to deal with all variations (eg, Lobatto, Radau and Legendre). LIBPFASST uses the augmented/virtual node route, and this PR goes with the square Q route. My guess is that, in the end, the logic for the square Q route will be more explicit and easier to follow. Let's see how it goes...

memmett · 2014-09-08T21:06:52Z

include/pfasst/encap/imex_sweeper.hpp

-          if (initial) { this->evaluate(0); }
+          if (initial) {
+            this->f_expl_eval(this->fs_expl_start, this->u_start, t + dt * this->quad->get_nodes()[0]);
+            this->f_impl_eval(this->fs_impl_start, this->u_start, t + dt * this->quad->get_nodes()[0]);


I don't think you want the t + dt... (just t). Also, I don't think the implicit piece is necessary here.

Somehow, I managed to "tune" the IMEX sweeper in such a way, that the implicit evaluation is necessary for overall convergence of the AdVec vanilla_sdc example. 😕

memmett · 2014-09-08T21:15:22Z

I think we should put this into a branch where we can all hack on it...

torbjoernk · 2014-09-09T06:39:06Z

I agree on having a common branch for that. I just pushed my current feature branch to upstream. You should all be able to work there. In case there are only us two working on that, it should be fine to directly push onto upstream without PRs.

I'll close this PR for now.

torbjoernk added 6 commits September 4, 2014 08:47

quad: the new and shiny quadrature

573de6c

quad: adjust tests

efbedb9

reordering quadrature

8011bd4

splitting into separate files in subdirectory note: just include 'quadrature.hpp' for all the quadrature stuff

tests: moving Polynomial tests into own file

3ce8f44

quad: a little template cleanup

606e291

tests: examples: first adjustments [skip ci]

5142110

torbjoernk added discussion labels Sep 4, 2014

torbjoernk added this to the Release 0.3.0 milestone Sep 4, 2014

torbjoernk added 2 commits September 4, 2014 10:50

quad: node-to-node deltas now member of quadrature

a6f92d5

quad: first adjustments to base sweepers

98230b2

torbjoernk force-pushed the feature/new-quadrature branch from 98230b2 to a6f92d5 Compare September 4, 2014 08:52

torbjoernk added 6 commits September 5, 2014 10:42

quad: remove unused parameter

1151354

quad: rename member variables and getters

8b61076

using naming convention of rest of codebase

quad: removing traits; minor fixes

3116d48

fixing GCC compilation

adc1f46

why does clang not complain? o.O

quad: removing traits & boilerplate code

9662023

everything should still work as before

sweeper: adjusting IMEX sweeper to new quadrature

6b826db

also adjusted AdVec Vanilla SDC example note: only vanilla_sdc compiles and runs

memmett reviewed Sep 5, 2014
View reviewed changes

sweeper: removing three node-methods and reinline the logic into sweep

5e6ccee

memmett reviewed Sep 8, 2014
View reviewed changes

torbjoernk added the delayed / obsolete label Sep 9, 2014

torbjoernk deleted the feature/new-quadrature branch September 11, 2014 07:07

This pull request was closed.

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

[WIP] Quadrature Revamping #86

[WIP] Quadrature Revamping #86

torbjoernk commented Sep 4, 2014

memmett commented Sep 4, 2014

pancetta commented Sep 4, 2014

mlminion commented Sep 4, 2014

pancetta commented Sep 4, 2014

memmett commented Sep 4, 2014

pancetta commented Sep 4, 2014

memmett commented Sep 4, 2014

mlminion commented Sep 4, 2014

pancetta commented Sep 4, 2014

mlminion commented Sep 4, 2014

torbjoernk commented Sep 5, 2014

memmett Sep 5, 2014

memmett commented Sep 5, 2014

memmett Sep 8, 2014

torbjoernk Sep 9, 2014

memmett commented Sep 8, 2014

torbjoernk commented Sep 9, 2014

[WIP] Quadrature Revamping #86

[WIP] Quadrature Revamping #86

Conversation

torbjoernk commented Sep 4, 2014

memmett commented Sep 4, 2014

pancetta commented Sep 4, 2014

mlminion commented Sep 4, 2014

pancetta commented Sep 4, 2014

memmett commented Sep 4, 2014

pancetta commented Sep 4, 2014

memmett commented Sep 4, 2014

mlminion commented Sep 4, 2014

pancetta commented Sep 4, 2014

mlminion commented Sep 4, 2014

torbjoernk commented Sep 5, 2014

memmett Sep 5, 2014

Choose a reason for hiding this comment

memmett commented Sep 5, 2014

memmett Sep 8, 2014

Choose a reason for hiding this comment

torbjoernk Sep 9, 2014

Choose a reason for hiding this comment

memmett commented Sep 8, 2014

torbjoernk commented Sep 9, 2014