TL: rework on node formulation tutorial

docs/conf.py

@@ -56,6 +56,7 @@
     "amsmath",
     "dollarmath",
 ]
+myst_heading_anchors = 2
 
 favicons = [
     "android-chrome-192x192.png",

docs/notebooks/01_qCoeffs.ipynb

@@ -349,14 +349,8 @@
   }
  ],
  "metadata": {
-  "kernelspec": {
-   "display_name": "env",
-   "language": "python",
-   "name": "python3"
-  },
   "language_info": {
-   "name": "python",
-   "version": "3.11.9"
+   "name": "python"
   }
  },
  "nbformat": 4,

docs/notebooks/11_nodeFormulation.ipynb

@@ -4,19 +4,17 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "# Extended Step 1 : Zero-to-Nodes (Z2N) and Node-to-Node (N2N) formulations\n",
+    "# Extended Step 1 : Zero-to-Nodes (Z2N) and Node-to-Node (N2N)\n",
     "\n",
-    "📜 _If you already know about SDC from the [original paper](https://link.springer.com/content/pdf/10.1023/A:1022338906936.pdf) of Dutt, Greengard and Rokhlin, you may notice that their description is very different from the one given [in Step 4](./04_sdc.ipynb) ..._\n",
+    "📜 _If you already know about SDC from the [original paper](https://link.springer.com/content/pdf/10.1023/A:1022338906936.pdf) of Dutt, Greengard & Rokhlin, you may notice that their description is very different from the one given [in Step 4](./04_sdc.ipynb) ..._\n",
     "\n",
     "Indeed, this tutorial introduced SDC using a **Zero-to-Nodes formulation (Z2N)**, which describes the SDC node updates from the initial step solution (zero) to the node. This approach is identical as the one used to describe Runge-Kutta methods with Butcher tables in the literature.\n",
     "\n",
     "The SDC authors however used a different formulation, namely the **Node-to-Node formulation (N2N)**, which describes the node update from one node to the next. While both formulations can produce identical SDC schemes, they have some fundamental differences from an implementation perspective, and leads to different generalizations of SDC.\n",
     "\n",
-    "We describe here how to switch from Z2N to N2N (and vice-versa), and how to implement the N2N formulation using `qmat`.\n",
+    "## Deriving the N2N formulation\n",
     "\n",
-    "## From Zero-to-Nodes to Node-to-Node\n",
-    "\n",
-    "Starting from the Z2N update on the Dahlquist problem :\n",
+    "We start from the Z2N update on the Dahlquist problem :\n",
     "\n",
     "$$\n",
     "u^{k+1} - \\lambda\\Delta{t}Q_\\Delta u^{k+1} = u_n + \\lambda\\Delta{t}(Q-Q_\\Delta)u^k,\n",
@@ -43,28 +41,32 @@
     "$$\n",
     "\n",
     "Now subtracting the update formula for $u^{k+1}_m$ from $u^{k+1}_{m+1}$,\n",
-    "we get the **generic N2N sweep formula** for $m > 0$ (starting from the second node) :\n",
+    "we get for $m > 0$ (starting from the second node) :\n",
     "\n",
     "$$\n",
     "u^{k+1}_{m+1} = u^{k+1}_m\n",
     "    + \\lambda\\Delta{t}\\sum_{j=1}^{m+1}\\left(q^\\Delta_{m+1,j} - q^\\Delta_{m,j}\\right)\\left(u^{k+1}_{j} - u^{k}_{j}\\right)\n",
     "    + \\lambda\\Delta{t}\\sum_{j=1}^{M}\\left(q_{m+1,j}-q_{m,j}\\right)u^{k}_{j},\n",
     "$$\n",
     "\n",
-    "and for the first node (no subtraction):\n",
+    "and for the first node (no difference) :\n",
     "\n",
     "$$\n",
     "u^{k+1}_{1} = u_n \n",
     "    + \\lambda\\Delta{t}q^\\Delta_{1,1} (u^{k+1}_{1} - u^{k}_{1}) \n",
-    "    + \\lambda\\Delta{t}\\sum_{j=1}^{M}q_{1,j}u^{k}_{j}.\n",
+    "    + \\lambda\\Delta{t}\\sum_{j=1}^{M}q_{1,j}u^{k}_{j}\n",
     "$$\n",
     "\n",
-    "> :bell: Note that $q^\\Delta_{m,m+1}=0$ because of the lower triangular nature of $Q_\\Delta$, so we can add this coefficient\n",
-    "> in the generic N2N sweep formula to simplify it.\n",
+    "We call this the **Node-to-node (N2N) update**, since each node update for $u^{k+1}_{m+1}$\n",
+    "explicitly depends on the previously computed update $u^{k+1}_{m}$ (and the correction terms $\\sum ...$),\n",
+    "except for the first node where it depends on the initial step solution $u_n$.\n",
+    "\n",
+    "> 🔔 Note that $q^\\Delta_{m,m+1}=0$ because of the lower triangular nature of $Q_\\Delta$, \n",
+    "> so we can add this coefficient in the generic N2N sweep formula to simplify it.\n",
     "\n",
     "Defining $s_{m+1,j} = q_{m+1,j}-q_{m,j} \\; \\forall m \\in \\{1, M-1\\}$ and $s_{1,j} = q_{1,j}$,\n",
     "we note $S$ the matrix built with the $(s)_{i,j}$ coefficients,\n",
-    "and write the **generic N2N formula** into matrix formulation :\n",
+    "and write the **generic N2N sweep update** into matrix formulation :\n",
     "\n",
     "$$\n",
     "\\begin{pmatrix}\n",
@@ -80,18 +82,72 @@
     "+ \\lambda\\Delta{t}S_\\Delta (u^{k+1}-u^{k}) + \\lambda\\Delta{t}S u^k,\n",
     "$$\n",
     "\n",
-    "where $S_\\Delta$ is built from the $Q_\\Delta$ matrix the same way as $S$ from $Q$.\n",
+    "where $S_\\Delta$ is built from the $Q_\\Delta$ matrix the same way as $S$ from $Q$. \n",
+    "Now you may remark that the N2N formula is actually more complex than the Z2N formula given above.\n",
+    "Thing is : some specific $Q_\\Delta$ coefficients allows to simplify it a lot ..."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Backward-Euler based sweep\n",
+    "\n",
+    "Consider now one of the original SDC form proposed by Dutt, Greengard & Rokhlin : the Backward-Euler based sweep.\n",
+    "The associated $Q_\\Delta$ coefficients are implemented in `qmat`, and considering a simple illustrative node distribution we obtain (using the default N2N formulation) :"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([[0.1, 0. , 0. , 0. ],\n",
+       "       [0.1, 0.2, 0. , 0. ],\n",
+       "       [0.1, 0.2, 0.4, 0. ],\n",
+       "       [0.1, 0.2, 0.4, 0.3]])"
+      ]
+     },
+     "execution_count": 1,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "from qmat import genQDeltaCoeffs\n",
+    "genQDeltaCoeffs(\"BE\", nodes=[0.1, 0.3, 0.7, 1.0])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "This is a triangular matrix with non-zero diagonal (implicit sweep), but with **all non-zero coefficients in each columns being identical**. This implies that for $m \\in \\{1, M\\}$ :\n",
+    "\n",
+    "$$\n",
+    "\\lambda\\Delta{t}\\sum_{j=1}^{m} s^\\Delta_{m+1,j}\\left(u^{k+1}_{j} - u^{k}_{j}\\right) = \n",
+    "\\lambda\\Delta{t}\\sum_{j=1}^{m}\\left(q^\\Delta_{m+1,j} - q^\\Delta_{m,j}\\right)\\left(u^{k+1}_{j} - u^{k}_{j}\\right) = 0\n",
+    "$$\n",
     "\n",
+    "thus simplifies the N2N formula into :\n",
     "\n",
-    "**Backward-Euler based sweep**\n",
+    "$$\n",
+    "u^{k+1}_{m+1} = u^{k+1}_m\n",
+    "    + \\lambda\\Delta{t}s^\\Delta_{m+1,m+1}\\left(u^{k+1}_{m+1} - u^{k}_{m+1}\\right)\n",
+    "    + \\lambda\\Delta{t}\\sum_{j=1}^{M}s_{m+1,j}u^{k}_{j},\n",
+    "$$\n",
     "\n",
-    "Consider now one of the original SDC form proposed by Dutt, Greengard and Rokhlin : the Backward-Euler based sweep.\n",
-    "The associated $Q_\\Delta$ coefficients are implemented in `qmat`, and considering a simple illustrative node distribution we obtain :"
+    "where we note $\\forall k\\;u^{k}_0 := u_n$, such that the formula can be applied on all nodes.\n",
+    "In fact, the $S_\\Delta$ matrix is diagonal, as we can see using `qmat` setting `form=\"N2N\"` for the \n",
+    "`genQDeltaCoeffs` function :"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": 2,
    "metadata": {},
    "outputs": [
     {
@@ -103,38 +159,167 @@
        "       [0. , 0. , 0. , 0.3]])"
       ]
      },
-     "execution_count": 1,
+     "execution_count": 2,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
-    "from qmat import genQDeltaCoeffs\n",
     "genQDeltaCoeffs(\"BE\", form=\"N2N\", nodes=[0.1, 0.3, 0.7, 1.0])"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "We obtain a triangular matrix with non-zero diagonal (implicit sweep), but with **all non-zero coefficients in each columns being identical**. This implies that for $m \\in \\{1, M\\}$ :\n",
+    "In this case in particular, we have $s^\\Delta_{m+1,m+1} = q^\\Delta_{m+1,m+1} = \\tau_{m+1}-\\tau_{m}$,\n",
+    "where $(\\tau_m)_{1\\leq m \\leq M} \\in [0, 1]$ are the nodes positions and $\\tau_0=0$.\n",
+    "This is then **exactly** the SDC formula for Backward Euler introduced by\n",
+    "[Dutt, Greengard & Rokhlin](https://link.springer.com/content/pdf/10.1023/A:1022338906936.pdf).\n",
+    "\n",
+    "> 🔔 Note that the **correction term** $\\lambda\\Delta{t}s^\\Delta_{m+1,m+1}\\left(u^{k+1}_{m+1} - u^{k}_{m+1}\\right)$\n",
+    "> depends on the result of the sweep update $u^{k+1}_{m+1}$, requiring then an **implicit** solver. "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Forward Euler based sweep\n",
+    "\n",
+    "An other sweep update was proposed by Dutt, Greengard & Rokhlin in the same spirit as for Backward Euler, \n",
+    "but this time using a Forward Euler based correction.\n",
+    "\n",
+    "Consider the $Q_\\Delta$ coefficients implemented in `qmat` for the same node distribution as before :"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([[0. , 0. , 0. , 0. ],\n",
+       "       [0.2, 0. , 0. , 0. ],\n",
+       "       [0.2, 0.4, 0. , 0. ],\n",
+       "       [0.2, 0.4, 0.3, 0. ]])"
+      ]
+     },
+     "execution_count": 3,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "genQDeltaCoeffs(\"FE\", nodes=[0.1, 0.3, 0.7, 1.0])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Those are similar to the ones obtained for Backward Euler, except they are \"shifted\" to the left and the diagonal contains only zeros. Looking at the corresponding $S_\\Delta$ coefficients, we obtain :\n",
     "\n",
     "$$\n",
-    "\\lambda\\Delta{t}\\sum_{i=1}^{m}\\left(q^\\Delta_{m+1,i} - q^\\Delta_{m,i}\\right)\\left(u^{k+1}_{i} - u^{k}_{i}\\right) = 0\n",
+    "u^{k+1}_{m+1} = u^{k+1}_m\n",
+    "    + \\lambda\\Delta{t}s^\\Delta_{m,m}\\left(u^{k+1}_{m} - u^{k}_{m}\\right)\n",
+    "    + \\lambda\\Delta{t}\\sum_{j=1}^{M}s_{m+1,j}u^{k}_{j},\n",
     "$$\n",
     "\n",
-    "which simplifies the N2N formula into :\n",
+    "with $s^\\Delta_{m,m}=q^\\Delta_{m,m}=\\tau_{m+1}-\\tau_{m}$ where $(\\tau_m)_{1\\leq m \\leq M} \\in [0, 1]$ \n",
+    "are the nodes positions and $\\tau_0=0$ (as before),\n",
+    "as we can see looking at the $S_\\Delta$ coefficient generated by `qmat` :"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([[0. , 0. , 0. , 0. ],\n",
+       "       [0.2, 0. , 0. , 0. ],\n",
+       "       [0. , 0.4, 0. , 0. ],\n",
+       "       [0. , 0. , 0.3, 0. ]])"
+      ]
+     },
+     "execution_count": 4,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "genQDeltaCoeffs(\"FE\", form=\"N2N\", nodes=[0.1, 0.3, 0.7, 1.0])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "This is similar as for Backward Euler, except the correction term does not depend\n",
+    "on $u^{k+1}_{m+1}$ this time, hence making this update fully **explicit**."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Additional notes\n",
+    "\n",
+    "Other type of sweep have a simplified form in N2N formulation, e.g using the trapezoid rule (Crank-Nicholson) : "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([[0.05, 0.  , 0.  , 0.  ],\n",
+       "       [0.1 , 0.1 , 0.  , 0.  ],\n",
+       "       [0.  , 0.2 , 0.2 , 0.  ],\n",
+       "       [0.  , 0.  , 0.15, 0.15]])"
+      ]
+     },
+     "execution_count": 5,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "genQDeltaCoeffs(\"TRAP\", form=\"N2N\", nodes=[0.1, 0.3, 0.7, 1.0])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "which produce a sweep update of this form :\n",
     "\n",
     "$$\n",
     "u^{k+1}_{m+1} = u^{k+1}_m\n",
-    "    + \\lambda\\Delta{t}q^\\Delta_{m+1,m+1}\\left(u^{k+1}_{m+1} - u^{k}_{m+1}\\right)\n",
-    "    + \\lambda\\Delta{t}\\sum_{j=1}^{M}\\left(q_{m+1,j}-q_{m,j}\\right)u^{k}_{j},\n",
+    "    + \\lambda\\Delta{t}s^\\Delta_{m+1,m+1}\\left(u^{k+1}_{m+1} - u^{k}_{m+1}\\right)\n",
+    "    + \\lambda\\Delta{t}s^\\Delta_{m,m}\\left(u^{k+1}_{m} - u^{k}_{m}\\right)\n",
+    "    + \\lambda\\Delta{t}\\sum_{j=1}^{M}s_{m+1,j}u^{k}_{j}.\n",
     "$$\n",
     "\n",
-    "where we can note $u^{k}_0 := u_n\\quad\\forall k$ so that the formula can be applied on all nodes.\n",
-    "Defining $s_{m+1,j} = q_{m+1,j}-q_{m,j} \\quad \\forall m \\in \\{1, M-1\\}$ and $s_{1,j} = q_{1,j}$,\n",
-    "we note $S$ the matrix built with the $(s)_{i,j}$ coefficients,\n",
-    "and \n"
+    "From an algorithmic perspective, implementing SDC into N2N form for Backward / Forward Euler or the trapezoid rule\n",
+    "is usually more efficient than the Z2N form, as it usually requires less floating point operations during sweep \n",
+    "(correction terms use all node solution in Z2N). However, the N2N formulation has two major inconvenient when\n",
+    "considering generic SDC methods :\n",
+    "\n",
+    "1. Only a few type of $Q_\\Delta$ coefficients allows a simplified N2N formulation, while some other (like LU for instance), don't have a simplified N2N formulation, hence making the Z2N implementation more efficient\n",
+    "2. The N2N formulation is inherently sequential, and do not allow the design of parallel SDC variant that have a diagonal $Q_\\Delta$ matrix, allowing to perform all node update in parallel across the nodes\n",
+    "\n",
+    "Furthermore, while the N2N formulation was historically used to describe SDC, many later SDC-related publications\n",
+    "have been using the Z2N formulation since it's more convenient for the analysis, when looking to SDC as a preconditioned\n",
+    "fixed point iteration."
    ]
   }
  ],