Subcycling.md

+# A Short Film About Subcycling
+
+This page contains a short description of how we implemented subcycling in the `cemracs14`
+branch. The algorithm is heavily inspired by the pseudo-code given by A.M. Khokhlov in
+the [Fully Threaded Tree](http://cds.cern.ch/record/318999/files/9701194.pdf) article. The
+algorithm is however slightly modified to take into account ghost cells between processes.
+
+The pseudo-code is:
+
+```
+#!python
+import antigravity
+
+# iterators.c:update_cell_subcycle_iter_volume_fn:335
+def update(l):
+    for all cells of level (l):
+        for all directions (x, y, z):
+            k = neighbor
+
+            if k is a leaf of level <= l:
+                compute flux at interface
+                update cell value
+
+    for all cells of level (l - 1):
+        for all directions (x, y, z):
+            k = neighbor
+
+            if k is a leaf of level l:
+                compute flux at interface
+                update cell value
+
+# solver.c:solver_subcycle:259
+def subcycle(level):
+    if level > MAX_LEVEL:
+        return
+
+    for all cells of level (level):
+        copy old value into new
+
+    subcycle(level + 1)
+    subcycle(level + 1)
+
+    dtlevel = dt * pow(2, MAX_LEVEL - level)
+
+    for all ghosts of level (level) and (level - 1):
+        update ghost values from neighboring processes
+
+    update(level)
+
+    for all cells of level (level):
+        copy new value into old
+```
+
+A few mentions are in order probably.
+
+`dt` is the time step as computed at the most refined level. To 
+get a time step at other, coarser, level we have supposed that
+the CFL condition is tightly tied to the mesh structure and
+since we known that `dx` gets divided by 2 when we increase in level,
+we have done the same thing for `dt`. In short, this means that
+the time step at a level `l < lmax` is:
+
+```
+dt(l) = dt * 2^(lmax - l)
+```
+
+The cell update routine is implemented in such a way that it does
+not allow computing the same flux twice. First we have a loop
+over cells of level l and we compute the flux with cells that are
+smaller or equal in level, ignoring the ones that are bigger.
+
+Secondly, we update cells of level (l - 1) with the flux coming from
+cells of level l. This flux should have been added to these cells when
+computing it from the point of view of the cell at level l in the
+previous loop.
+
+The second loop is necessary to correctly deal with ghost cells between
+processes (probably not the most efficient solution). The problem
+it tries to solve is when we have a cell of level (l - 1) which has
+one or two neighbors of level l that are ghost cells. This cell would
+not be update when cycling over cells of level l, because we only loop
+over real cells, so it it necessary to also loop over cells of level
+(l - 1) to cover this edge case.
\ No newline at end of file