Skip to content

Commit

Permalink
minor corrections
Browse files Browse the repository at this point in the history
  • Loading branch information
Ig-dolci committed Feb 9, 2024
1 parent 21bec1c commit eb08589
Showing 1 changed file with 24 additions and 24 deletions.
48 changes: 24 additions & 24 deletions docs/notebooks/burger.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -57,24 +57,25 @@
"where $v \\in V$ is an arbitrary test function. The Jacobian of $F(u^{n+1}, u^n, v)$ is obtained through the Gateaux derivative. \n",
"\n",
"The time sequence of the adjoint system for the functional (Eq. (1)) is given by:\n",
"\n",
"\\begin{bmatrix}\n",
" \\frac{\\partial F_{N + 1}}{\\partial u^{N + 1}} \\lambda^{N + 1} + \\frac{\\partial I}{\\partial u^{N + 1}} = 0 \\\\\n",
"$$\n",
"\\begin{align*}\n",
" \\frac{\\partial F_{N + 1}}{\\partial u^{N + 1}} \\lambda^{N + 1} &= \\frac{\\partial I}{\\partial u^{N + 1}} \\\\\n",
" \\\\\n",
" \\frac{\\partial F_{N}}{\\partial u^{N}} \\lambda^{N} + \\frac{\\partial F_{N + 1}}{\\partial u^{N}} \\lambda^{N + 1} = 0 \\\\\n",
" \\frac{\\partial F_{N}}{\\partial u^{N}} \\lambda^{N} + \\frac{\\partial F_{N + 1}}{\\partial u^{N}} \\lambda^{N + 1} &= 0 \\\\\n",
" \\\\\n",
" \\vdots \\\\\n",
" \\\\\n",
" \\frac{\\partial F_{1}}{\\partial u^{1}} \\lambda^{1} + \\frac{\\partial F_{2}}{\\partial u^{1}} \\lambda^{2} = 0 \\\\\n",
" \\frac{\\partial F_{1}}{\\partial u^{1}} \\lambda^{1} + \\frac{\\partial F_{2}}{\\partial u^{1}} \\lambda^{2} &= 0 \\\\\n",
" \\\\\n",
" \\frac{\\partial F_{0}}{\\partial u^{0}} \\lambda^{0} + \\frac{\\partial F_{1}}{\\partial u^{0}} \\lambda^{1} = 0\n",
"\\end{bmatrix} \\tag{5}\n",
"\n",
" \\frac{\\partial F_{0}}{\\partial u^{0}} \\lambda^{0} + \\frac{\\partial F_{1}}{\\partial u^{0}} \\lambda^{1} &= 0\n",
"\\end{align*}\n",
"\\tag{5}\n",
"$$\n",
"where $\\lambda^{n}$ is the adjoint variable at time step $n$.\n",
"\n",
"The adjoint system is solved backward in time, we can verify this by taking the adjoint of Eq. (5), where at the time step zero we have two unknowns, $\\lambda^{0}$ and $\\lambda^{1}$, and at the time step $N$ we have only one unknown, $\\lambda^{N}$.\n",
"\n",
"For additonal details on the adjoint-based gradient computation, you can refer to the following references [1, 2].\n"
"For additional details on the discrete adjoint-based gradient computation, you may refer to the following references [1, 2].\n"
]
},
{
Expand All @@ -87,18 +88,18 @@
"The Burger's equation is discretised using the Finite Element Method (FEM). We use continuous first-order Lagrange basis functions to discretise the spatial domain.\n",
"\n",
"After discretising, the residual of the forward system is given by:\n",
"$$F = M \\mathbf{U}^{n + 1} - M \\mathbf{U}^{n} + \\Delta t \\left(C(\\mathbf{U}^{n + 1}) \\mathbf{U}^{n + 1} - K \\mathbf{U}^{n + 1}\\right) = 0, \\tag{6}$$\n",
"$$\\mathbf{F} = \\mathbf{M} \\mathbf{U}^{n + 1} - \\mathbf{M} \\mathbf{U}^{n} + \\Delta t \\left(\\mathbf{C}(\\mathbf{U}^{n + 1}) \\mathbf{U}^{n + 1} - \\mathbf{K} \\mathbf{U}^{n + 1}\\right) = 0, \\tag{6}$$\n",
"where $M$ and $K$ are the mass and stiffness, respectively. $U$ represents the solution vector. The matrices $C$ is obtained from the advection term. The superscript $n$ denotes the time step, and $\\Delta t$ is the time step size. The solution vector $\\mathbf{U}^{n + 1}$ is obtained by solving the non-linear forward system using the Newton method. The Jacobian of the forward system has the form:\n",
"\n",
"Let us consider the residual of the forward system as $R(\\mathbf{U}^{n + 1}, \\mathbf{U}^{n})$. The Jacobian of the forward system is given by:\n",
"$$J_{n, n} = \\frac{\\partial F^{n}}{\\partial U^{n}} = M + \\Delta t \\left(J_C(\\mathbf{U}^{n}) - K \\right), \\tag{7}$$\n",
"$$J_{n + 1, n} = \\frac{\\partial F^{n + 1}}{\\partial U^{n}} = - M, \\tag{8}$$\n",
"where $J_C$ is the Jacobian of the advection term. Considering the adjoint as shownd in Eq. (5), and the Jacobian of the forward system (Eq. (7) a (8)), we arrive at the following expression for the adjoint system.\n",
"$$J^{T}_{n, n} \\Lambda^{n} + J^{T}_{n+1, n} \\Lambda^{n + 1} = 0, \\tag{9}$$\n",
"Let us consider the residual of the forward system as $\\mathbf{R}(\\mathbf{U}^{n + 1}, \\mathbf{U}^{n})$. The Jacobian of the forward system is given by:\n",
"$$\\mathbf{J}_{n, n} = \\frac{\\partial \\mathbf{F}^{n}}{\\partial U^{n}} = \\mathbf{M} + \\Delta t \\left(\\mathbf{J_C}(\\mathbf{U}^{n}) - \\mathbf{K} \\right), \\tag{7}$$\n",
"$$\\mathbf{J_{n + 1, n}} = \\frac{\\partial \\mathbf{F}^{n + 1}}{\\partial \\mathbf{U}^{n}} = - \\mathbf{M}, \\tag{8}$$\n",
"where $\\mathbf{J_C}$ is the Jacobian of the advection term. Considering the adjoint as shownd in Eq. (5), and the Jacobian of the forward system (Eq. (7) a (8)), we arrive at the following expression for the adjoint system.\n",
"$$\\mathbf{J}^{T}_{n, n} \\Lambda^{n} + \\mathbf{J}^{T}_{n+1, n} \\boldsymbol{\\Lambda}^{n + 1} = 0, \\tag{9}$$\n",
"for $n = N - 1, \\ldots, 0$. $\\Lambda$ is the adjoint solution vector. \n",
"\n",
"At the time $n = N$, the adjoint is initialized by solving the following system:\n",
"$$J^{T}_{N,N} \\Lambda^{N} + \\frac{\\partial I}{\\partial U^{N}} = 0. \\tag{10}$$\n"
"$$\\mathbf{J}^{T}_{N,N} \\boldsymbol{\\Lambda}^{N} + \\frac{\\partial \\mathbf{I}}{\\partial \\mathbf{U}^{N}} = 0. \\tag{10}$$\n"
]
},
{
Expand Down Expand Up @@ -171,7 +172,7 @@
" if write_adj_deps:\n",
" self._store_data(u, step, storage, write_adj_deps, write_ics)\n",
" def _residual(u_new):\n",
" C = self._convection_matrix(u_new)\n",
" C = self._advection_matrix(u_new)\n",
" residual = (M + self.model[\"dt\"] * (- K + C)).dot(u_new) - M.dot(u)\n",
" residual[0] = residual[self.model[\"nx\"] - 1] = 0\n",
" return residual\n",
Expand Down Expand Up @@ -218,12 +219,6 @@
" step -= 1\n",
" self.adjoint_work_memory[StorageType.WORK][step] = copy.deepcopy(u)\n",
"\n",
" def _initialize_adjoint(self):\n",
" J = self._jacobian(self.forward_final_solution, self.model[\"max_n\"], self.model[\"max_n\"])\n",
" M = self._mass_matrix()\n",
" self.adjoint_work_memory[StorageType.WORK][self.model[\"max_n\"]] = spsolve(J.T, - self.forward_final_solution)\n",
"\n",
"\n",
" def gradient(self):\n",
" \"\"\"Compute the adjoint-based gradient.\n",
"\n",
Expand Down Expand Up @@ -276,6 +271,11 @@
" if step in self.restart_forward[from_storage]:\n",
" del self.restart_forward[from_storage][step]\n",
"\n",
" def _initialize_adjoint(self):\n",
" J = self._jacobian(self.forward_final_solution, self.model[\"max_n\"], self.model[\"max_n\"])\n",
" M = self._mass_matrix()\n",
" self.adjoint_work_memory[StorageType.WORK][self.model[\"max_n\"]] = spsolve(J.T, - self.forward_final_solution)\n",
"\n",
" def _jacobian(self, u, i, j):\n",
" M = self._mass_matrix()\n",
" if i == j:\n",
Expand Down Expand Up @@ -327,7 +327,7 @@
" K[i, i] = - 2 * b\n",
" return K\n",
" \n",
" def _convection_matrix(self, u):\n",
" def _advection_matrix(self, u):\n",
" num_nodes = self.model[\"nx\"]\n",
" # 1D mesh is uniform. Thus, the mesh spacing is constant.\n",
" h = self.model[\"lx\"] / self.model[\"nx\"] \n",
Expand Down

0 comments on commit eb08589

Please sign in to comment.