SDC Block Jacobi

[JensHahne]

sdc_block_jacobi = BlockIteration(
    "phid*h*u_{n}^k + (i-phid*phi)* u_{n+1}^{k}",
    "i * u_{n}^k",
    i=one, phid = BlockOperator('$\Phi_{\delta}^{-1}$'), phi = BlockOperator('$\Phi$'), h = BlockOperator('$\chi$'))

PintRun(pararealSC_new, nBlocks, None)

[tlunet]

I've debbuged the copy method of the BlockOperator class, that was missing one exception for the prediction rule.

However, there is now an other error on the graph side, see https://github.com/Parallel-in-Time/time4apint/blob/gaia/notebook/02_Elementary_primary_block_iterations.ipynb

[JensHahne]

I fixed the bug, somehow it was not clear from the graph side how to handle the identity. I have updated the notebook with two versions:

1. The first version is the original version. But the problem is that the optimization works maybe a bit too well here and optimizes everything out of the bottom right of the graph. From my side the graph looks unexpected but somehow correct.
2. Instead of the operator one for the identical matrix, I used an operator I. In this case the lower right part is not optimized. This looks a bit ugly for the prediction phase, but actually everything is correct. Since we minimize the tasks, we use for example the already existing u_2^0=I * I *u_0^0 for computing u_1^1 (which also requires the operation I * I * u_0^0).

In both cases there is a "problem". From the PFASST point of view, each time block n would be located on a different process, so we don't need communication for all iterations of the SDC block jacobi. But with the mathematically optimized version, we overlook the parallel aspect of the same cheap operation being performed on multiple processes. I don't know if this is really a problem in the end, but the optimization doesn't see these parallel ideas at this point. We should keep that in mind and think about it.

One idea might be to optimize only those tasks that have a real impact, i.e., where the cost is above a certain threshold. But this is just a spontaneous idea, I need to think more about it

[tlunet]

Actually, the original version is not really wrong : since you only use a copy for the initial guess, the sympy simplification does recognize that you do a copy of the initial solution and thus, the "real" algorithm is purely sequential (even if, you could have done some computation with the copied initial solution, but those would have been redundant). I believe this is equivalent to the simplification done for Parareal when $k=n$.

It appears that the predictor step is very important in the graph generation (which is not too surprising), hence the solution in figure 2 seems more general (even if the 2* factor is true because of how we built the predictor, but it still looks a bit weird). So in the context of fully analyzing the different tasks for the algorithm, I think it's still ok ... you should just know that block jacobi itself is a smoother, so not really used as standalone algorithm.

One improvement that could be good though : a way to plot the graph from non-dependent $u_{n}^0$, i.e without specifying a predictor step, so there would be as much starting nodes as there is blocks (that is, the $u_{n}^0$ for $n\in{1, ..., N_{blocks}}$). Could you implement something like this ? This could allow to analyze the block iteration independent from the predictor (for instance, if you use a randomized initial guess).

Also, could you represent the $u_n^0$ nodes for any cases (even when the predictor step is defined) ?
An put all the $n_n^k$ nodes at their corresponding edge of the grid (i.e $u_n^{k}$ at plot coordinate $(x,y)=(n,k)$ ?

[tlunet]

PS : changed the notebook name to ApproximateBlockJacobi, since we are specifically focusing on Approximate Block Jacobi here ...

I've added the notebook n°3 with all three elemental primary block iterations : 03_PrimaryBlockIteration. Eventually, if you implement the graph without predictor, this could be interesting to remove the predictor step, in order to compare primary block iterations only (in particular, there parallel efficiency considering equivalent block operator cost ...

[JensHahne]

One improvement that could be good though : a way to plot the graph from non-dependent $u_n^0$, i.e without specifying a predictor step, so there would be as much starting nodes as there is blocks (that is, the $u_n^0$ for $n \in 1,...,N_{blocks}$). Could you implement something like this ? This could allow to analyze the block iteration independent from the predictor (for instance, if you use a randomized initial guess).

The predictor is now optional. I have added examples without predictor to 03_PrimaryBlockIteration. But this is only handled on the graph side, maybe it needs to be adjusted on the BlockIteration side for convergence analysis (maybe not).

An put all the $n_n^k$ nodes at their corresponding edge of the grid (i.e $u_n^{k}$ at plot coordinate $(x,y)=(n,k)$ ?

Done.

Also, could you represent the u_n^0 nodes for any cases (even when the predictor step is defined) ?

I don't think that's a good idea. In the end, the node representation of a task is somehow not good at the moment. Each task represents both a result and an operation such as $u_n^{k+1}=\chi u_n^k$. Currently, we sometimes represent a node by its result and sometimes by (parts of) its operation. I could just say that each node $u_n^0$ is represented by its result, but I don't think that would improve readability, since we would then lose parts of the operation. Probably the best way would be to represent each task by result=operation, but that's just too long. I'll try to play around a bit with the sizes of each node/font size, maybe we'll find a way to represent a task by both.

Also, we should keep in mind that the graph representation is not really the final goal. The graph is very nice for small examples and to understand if everything looks as expected, but I think in the end we are more interested in the minimum possible runtime of a variant, which can be determined by the graph (so for larger applications you would probably never draw the graph, but only the minimum runtime).