Skip to content

Commit

Permalink
TL: update notebooks
Browse files Browse the repository at this point in the history
  • Loading branch information
tlunet committed Jul 19, 2024
1 parent e48429e commit 6d4c64e
Show file tree
Hide file tree
Showing 4 changed files with 153 additions and 63 deletions.
4 changes: 2 additions & 2 deletions docs/notebooks/02_rk.ipynb

Large diffs are not rendered by default.

4 changes: 2 additions & 2 deletions docs/notebooks/04_sdc.ipynb

Large diffs are not rendered by default.

148 changes: 148 additions & 0 deletions docs/notebooks/11_nodeFormulation.ipynb
Original file line number Diff line number Diff line change
@@ -0,0 +1,148 @@
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Extended Step 1 : Zero-to-Nodes (Z2N) and Node-to-Node (N2N) formulations\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",
"\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",
"\n",
"## From Zero-to-Nodes to Node-to-Node\n",
"\n",
"Starting 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",
"$$\n",
"\n",
"we can expend it to an update for each node solution $u_{m} \\simeq u(t_n+\\tau_m\\Delta{t}),\\; m \\in \\{1,\\dots,M\\}$:\n",
"\n",
"$$\n",
"u^{k+1}_{m+1} - \\lambda\\Delta{t}\\sum_{j=1}^{m+1}q^\\Delta_{m+1,j} u^{k+1}_{j} \n",
" = u_n \n",
" + \\lambda\\Delta{t}\\sum_{j=1}^{M}q_{m+1,j}u^{k}_{j}\n",
" - \\lambda\\Delta{t}\\sum_{j=1}^{m+1}q^\\Delta_{m+1,j}u^{k}_{j},\n",
"$$\n",
"\n",
"where $u_n$ is the initial solution for the time-step (scalar, abusing notation again ...),\n",
"and we note $(q^\\Delta)_{i,j} := Q_\\Delta$ and $(q)_{i,j} := Q$.\n",
"Rearranging and regrouping terms, we can write it like this :\n",
"\n",
"$$\n",
"u^{k+1}_{m+1} =\n",
" u_n \n",
" + \\lambda\\Delta{t}\\sum_{j=1}^{m+1}q^\\Delta_{m+1,j} (u^{k+1}_{j} - u^{k}_{j}) \n",
" + \\lambda\\Delta{t}\\sum_{j=1}^{M}q_{m+1,j}u^{k}_{j}.\n",
"$$\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",
"\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",
"\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",
"$$\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",
"\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",
"\n",
"$$\n",
"\\begin{pmatrix}\n",
"1 \\\\\n",
"-1 & 1 \\\\\n",
"& \\ddots & \\ddots \\\\\n",
"& & -1 & 1\n",
"\\end{pmatrix}\n",
"(u^{k+1} - u^{k}) \n",
"= \\begin{pmatrix}\n",
"u0 \\\\ 0 \\\\ \\vdots \\\\ 0 \n",
"\\end{pmatrix}\n",
"+ \\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",
"\n",
"\n",
"**Backward-Euler based sweep**\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 :"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[0.1, 0. , 0. , 0. ],\n",
" [0. , 0.2, 0. , 0. ],\n",
" [0. , 0. , 0.4, 0. ],\n",
" [0. , 0. , 0. , 0.3]])"
]
},
"execution_count": 1,
"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",
"\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",
"$$\n",
"\n",
"which simplifies the N2N formula into :\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",
"$$\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"
]
}
],
"metadata": {
"language_info": {
"name": "python"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
60 changes: 1 addition & 59 deletions docs/notebooks/11_prolongation.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -6,65 +6,7 @@
"source": [
"# Step 5 : generalizing the prolongation for RK-type and SDC-type time-steppers\n",
"\n",
"🛠️ In construction ...\n",
"\n",
"## Additional coefficients from $G$-generators\n",
"\n",
"While the `genQCoeffs` function and `genCoeffs` method of the $Q$-generator provide per default only the nodes, weights and $Q$ matrix (zero-to-nodes), \n",
"you can also use them to retrieve :\n",
"\n",
"- the $S$ matrix (node-to-node)\n",
"- the interpolation coefficients for the end-interval update `hCoeffs`"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"from qmat import genQCoeffs\n",
"from qmat.qcoeff.collocation import Collocation\n",
"coll = Collocation(nNodes=4, nodeType=\"LEGENDRE\", quadType=\"RADAU-RIGHT\")\n",
"\n",
"nodes, weights, Q, S = coll.genCoeffs(withS=True)\n",
"nodes, weights, Q, hCoeffs = genQCoeffs(\"coll\", hCoeffs=True, nNodes=4, nodeType=\"LEGENDRE\", quadType=\"RADAU-RIGHT\")\n",
"nodes, weights, Q, S, hCoeffs = coll.genCoeffs(withS=True, hCoeffs=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that $S$ is always returned before `hCoeffs`, if asked. Also, you can retrieve those directly from the generator object :"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"S :\n",
"[[ 0.08696371 -0.02660418 0.01262746 -0.00355515]\n",
" [ 0.10115441 0.18964047 -0.04050789 0.01029065]\n",
" [-0.0209262 0.19091672 0.19091672 -0.0209262 ]\n",
" [ 0.01029065 -0.04050789 0.18964047 0.10115441]]\n",
"hCoeffs : [-0.1139172 0.40076152 -0.81363245 1.52678813]\n"
]
}
],
"source": [
"coll = Collocation(nNodes=4, nodeType=\"LEGENDRE\", quadType=\"GAUSS\")\n",
"S = coll.S\n",
"hCoeffs = coll.hCoeffs\n",
"\n",
"print(\"S :\")\n",
"print(S)\n",
"print(\"hCoeffs :\", hCoeffs)"
"🛠️ In construction ..."
]
}
],
Expand Down

0 comments on commit 6d4c64e

Please sign in to comment.