diff --git a/course-rom/_attachments/homework/problem-set-1.ipynb b/course-rom/_attachments/homework-2023/problem-set-1.ipynb
similarity index 100%
rename from course-rom/_attachments/homework/problem-set-1.ipynb
rename to course-rom/_attachments/homework-2023/problem-set-1.ipynb
diff --git a/course-rom/_attachments/homework/problem-set-2.ipynb b/course-rom/_attachments/homework-2023/problem-set-2.ipynb
similarity index 100%
rename from course-rom/_attachments/homework/problem-set-2.ipynb
rename to course-rom/_attachments/homework-2023/problem-set-2.ipynb
diff --git a/course-rom/_attachments/homework/problem-set-3.ipynb b/course-rom/_attachments/homework-2023/problem-set-3.ipynb
similarity index 100%
rename from course-rom/_attachments/homework/problem-set-3.ipynb
rename to course-rom/_attachments/homework-2023/problem-set-3.ipynb
diff --git a/course-rom/_attachments/homework-2024-data.zip b/course-rom/_attachments/homework-2024-data.zip
new file mode 100644
index 0000000..6e5fbd8
Binary files /dev/null and b/course-rom/_attachments/homework-2024-data.zip differ
diff --git a/course-rom/_attachments/homework-2024/problem-set-1.ipynb b/course-rom/_attachments/homework-2024/problem-set-1.ipynb
new file mode 100644
index 0000000..24afc16
--- /dev/null
+++ b/course-rom/_attachments/homework-2024/problem-set-1.ipynb
@@ -0,0 +1 @@
+{"cells":[{"cell_type":"markdown","source":["# Problem Set 1: RB for Linear Affine Elliptic Problems\n\n","We consider the problem of designing a thermal fin to effectively remove heat from a surface. The two-dimensional fin, shown in [Figure 1](#fig:1), consists of a vertical central post and four horizontal subfins; the fin conducts heat from a prescribed uniform flux source at the root, $\\Gamma_{\\mathrm{root}}$ , through the large-surface-area subfins to surrounding flowing air. The fin is characterized by a five-component parameter vector, or input, latexmath:[\\mu_\n","= (\\mu_1 , \\mu_2, \\ldots, \\mu_5 )],where latexmath:[\\mu_i = k^i , i = 1, \\ldots\n",", 4], and $\\mu_5 = \\mathrm{Bi}$; $\\mu$ may take on any value in a specified design set $D \\subset \\mathbb{R}^5$.\n"],"metadata":{}},{"cell_type":"markdown","source":["\n","![image](fin.png)\n","\n"],"metadata":{"node_name":"image"}},{"cell_type":"markdown","source":["Here $k^i$ is the thermal conductivity of the ith subfin (normalized relative to the post conductivity $k^0 \\equiv 1$); and $\\mathrm{Bi}$ is the Biot number, a nondimensional heat transfer coefficient reflecting convective transport to the air at the fin surfaces (larger $\\mathrm{Bi}$ means better heat transfer). For example, suppose we choose a thermal fin with $k^1 = 0.4, k^2 = 0.6, k^3 = 0.8, k^4 = 1.2$, and $\\mathrm{Bi} = 0.1$; for this particular configuration $\\mu = \\{0.4, 0.6, 0.8, 1.2, 0.1\\}$, which corresponds to a single point in the set of all possible configurations D (the parameter or design set). The post is of width unity and height four; the subfins are of fixed thickness $t = 0.25$ and length $L = 2.5$.\n"],"metadata":{}},{"cell_type":"markdown","source":["We are interested in the design of this thermal fin, and we thus need to look at certain outputs or cost-functionals of the temperature as a function of $\\mu$. We choose for our output $T_{\\mathrm{root}}$, the average steady-state temperature of the fin root normalized by the prescribed heat flux into the fin root. The particular output chosen relates directly to the cooling efficiency of the fin — lower values of $T_{\\mathrm{root}}$ imply better thermal performance. The steady–state temperature distribution within the fin, $u(\\mu)$, is governed by the elliptic partial differential equation\n"],"metadata":{}},{"cell_type":"markdown","source":["\n$$\n\\label{eq:1}\n-k^i \\Delta u^i = 0 \\text{ in } \\Omega^i , i = 0, \\ldots, 4,\n$$\n"],"metadata":{}},{"cell_type":"markdown","source":["where $\\Delta$ is the Laplacian operator, and $u_i$ refers to the restriction of $u \\text{ to } \\Omega^i$ . Here $\\Omega^i$ is the region of the fin with conductivity $k^i , i = 0,\\ldots, 4: \\Omega^0$ is thus the central post, and $\\Omega^i , i = 1,\\ldots, 4$, corresponds to the four subfins. The entire fin domain is denoted $\\Omega (\\bar{\\Omega} = \\cup_{i=0}^4 \\bar{\\Omega}^i )$; the boundary $\\Omega$ is denoted $\\Gamma$. We must also ensure continuity of temperature and heat flux at the conductivity– discontinuity interfaces $\\Gamma^i_{int} \\equiv \\partial\\Omega^0 \\cap \\partial\\Omega^i , i = 1,\\ldots, 4$, where $\\partial\\Omega^i$ denotes the boundary of $\\Omega^i$, we have on $\\Gamma^i_{int} i = 1,\\ldots, 4$ :\n"],"metadata":{}},{"cell_type":"markdown","source":["\n$$\n\\begin{aligned}\n u^0 &= u^i \\\\\n -(\\nabla u^0 \\cdot n^i ) &= -k^i (\\nabla u^i \\cdot n^i )\n\\end{aligned}\n$$\n"],"metadata":{}},{"cell_type":"markdown","source":["here $n^i$ is the outward normal on $\\partial\\Omega^i$ . Finally, we introduce a Neumann flux boundary condition on the fin root\n"],"metadata":{}},{"cell_type":"markdown","source":["\n$$\n-(\\nabla u^0 \\cdot n^0 ) = -1 \\text{ on } \\Gamma_{\\mathrm{root}} ,\n$$\n"],"metadata":{}},{"cell_type":"markdown","source":["which models the heat source; and a Robin boundary condition\n"],"metadata":{}},{"cell_type":"markdown","source":["\n$$\n-k^i (\\nabla u^i \\cdot n^i ) = \\mathrm{Bi} u^i \\text{ on } \\Gamma^i_{ext} , i = 0,\\ldots, 4,\n$$\n"],"metadata":{}},{"cell_type":"markdown","source":["which models the convective heat losses. Here $\\Gamma^i_{ext}$ is that part of the boundary of $\\Omega^i$ exposed to the flowing fluid; note that $\\cup_{i=0}^4 \\Gamma^i_{ext} = \\Gamma\\backslash\\Gamma_{\\mathrm{root}}$. The average temperature at the root, $T_{\\mathrm{root}} (\\mu)$, can then be expressed as $\\ell^O(u(\\mu))$, where\n"],"metadata":{}},{"cell_type":"markdown","source":["\n$$\n\\ell^O (v) = \\int_{\\Gamma_{\\mathrm{root}}} v\n$$\n"],"metadata":{}},{"cell_type":"markdown","source":["(recall $\\Gamma_{\\mathrm{root}}$ is of length unity). Note that $\\ell(v) = \\ell^O(v)$ for this problem.\n\n","## Part 1 - Finite Element Approximation\n\n","We saw in class that the reduced basis approximation is based on a truth finite element approximation of the exact (or analytic) problem statement. To begin, we have to show that the exact problem described above does indeed satisfy the affine parameter dependence and thus fits into the framework shown in class.\n","a) Show that $u^e (\\mu) \\in X^e \\equiv H^1 (\\Omega)$ satisfies the weak form\n","\n$$\na(u^e (\\mu), v; \\mu) = \\ell(v), \\forall v \\in X^e ,\n$$\n","with\n","\n$$\n\\begin{aligned}\na(w, v; \\mu) &=\\sum_{i=0}^4 k^i \\int_{\\Omega^i} \\nabla w \\cdot \\nabla v dA\n+ \\mathrm{Bi}\\int_{\\Gamma\\backslash\\Gamma_{\\mathrm{root}}} w v dS,\\\\\n\\ell(v) &= \\int_{\\Gamma_{\\mathrm{root}}} v\n\\end{aligned}\n$$\n","b) Show that $u^e (\\mu)$ is the argument that minimizes\n","\n$$\n\\label{eq:2}\nJ(w) = \\frac{1}{2}\\sum_{i=0}^4 k^i \\int_{\\Omega_i} \\nabla w \\cdot\n\\nabla w dA + \\frac{\\mathrm{Bi}}{2}\n\\int_{\\Gamma\\backslash\\Gamma_{\\mathrm{root}}} w^2 dS -\n\\int_{\\Gamma_{\\mathrm{root}}} w dS\n$$\n","over all functions $w$ in $X^e$ .\n","We now consider the linear finite element space\n","\n$$\nX^{\\mathcal{N}} = \\{v \\in H^1 (\\Omega)| v|_{T_h} \\in \\mathbb{P}^1 (\\mathcal{T}_h ), \\forall T_h \\in \\mathcal{T}_h \\},\n$$\n","and look for $u^{\\mathcal{N}} (\\mu) \\in X^{\\mathcal{N}}$ such that\n","\n$$\na(u^{\\mathcal{N}} (\\mu), v; \\mu) = \\ell(v), \\forall v \\in X^{\\mathcal{N}} ;\n$$\n","our output of interest is then given by\n","\n$$\nT_{\\mathrm{root}}^\\mathcal{N}(\\mu) = \\ell^O (u^{\\mathcal{N}} (\\mu)).\n$$\n","Applying our usual nodal basis, we arrive at the matrix equations\n","\n$$\n\\begin{aligned}\n%\n A^{\\mathcal{N}} (\\mu) %\n u^{\\mathcal{N}} (\\mu) & = %\n F^{\\mathcal{N}} ,\\\\\nT_{\\mathrm{root}}^\\mathcal{N}(\\mu) &= (%\n L^{\\mathcal{N}} )^T %\n u^{\\mathcal{N}} (\\mu),\n\\end{aligned}\n$$\n","where latexmath:[%\n"," A^{\\mathcal{N}} \\in\n","\\mathbb{R}^{\\mathcal{N}\\times\\mathcal{N}} ,\n","%\n"," u^{\\mathcal{N}} \\in \\mathbb{R}^{\\mathcal{N}} ,\n","%\n"," F^{\\mathcal{N}} \\in \\mathbb{R}^{\\mathcal{N}} , and\n","%\n"," L^{\\mathcal{N}} \\in \\mathbb{R}^{\\mathcal{N}}]; here $\\mathcal{N}$ is the dimension of the finite element space $X^{\\mathcal{N}}$, which (given our natural boundary conditions) is equal to the number of nodes in $\\mathcal{T}_h$.\n\n","## Part 2 - Reduced-Basis Approximation\n\n","In general, the dimension of the finite element space, latexmath:[\\mathop{\\mathrm{dim}}X =\n","\\mathcal{N}], will be quite large (in particular if we were to treat the more realistic three-dimensional fin problem), and thus the solution of latexmath:[%\n"," A^{\\mathcal{N}} %\n"," u^{\\mathcal{N}} (\\mu) = %\n"," F^\\mathcal{N}] can be quite expensive. We thus investigate the reduced-basis methods that allow us to accurately and very rapidly predict $T_{\\mathrm{root}} (\\mu)$ in the limit of many evaluations — that is, at many different values of $\\mu$ — which is precisely the limit of interest in design and optimization studies. To derive the reduced-basis approximation we shall exploit the energy principle,\n","\n$$\nu(\\mu) = \\mathop{\\mathrm{argmin}}_{w \\in X} J(w),\n$$\n","where $J(w)$ is given by (#eq:2).\n","To begin, we introduce a sample in parameter space, latexmath:[S_N = \\{\\mu_1 ,\n"," \\mu_2 ,\\ldots, \\mu_N \\}] with $N \\ll \\mathcal{N}$. Each latexmath:[\\mu_i , i =\n","1,\\ldots, N] , belongs in the parameter set $\\mathcal{D}$. For our parameter set we choose $\\mathcal{D} = [0.1, 10.0$^4 \\times [0.01, 1.0]], that is, $0.1 \\leq k^i \\leq 10.0, i = 1,\\ldots, 4$ for the conductivities, and latexmath:[0.01\n","\\leq \\mathrm{Bi} \\leq 1.0] for the Biot number. We then introduce the reduced-basis space as\n","\n$$\nW_N = \\mathop{\\mathrm{span}}\\{u^\\mathcal{N} (\\mu_1 ), u^\\mathcal{N} (\\mu_2 ),\\ldots, u^\\mathcal{N}(\\mu_N ) \\}\n$$\n","where $u_N (\\mu_i )$ is the finite-element solution for latexmath:[\\mu =\n","\\mu_i].\n","To simplify the notation, we define $\\xi^i \\in X$ as latexmath:[\\xi^i = u_N\n","(\\mu_i ), i = 1,\\ldots, N]; we can then write latexmath:[W_N = span\\{\\xi^i , i =\n"," 1,\\ldots, N \\}]. Recall that $W_N = \\mathop{\\mathrm{span}}\\{\\xi^i , i = 1,\\ldots, N \\}$ means that $W_N$ consists of all functions in $X$ that can be expressed as a linear combination of the $\\xi^i$ ; that is, any member $v_N$ of $W_N$ can be represented as\n","\n$$\nv_N =\\sum_{i=1}^N \\beta^j \\xi^j ,\n$$\n","for some unique choice of latexmath:[\\beta^j \\in \\mathbb{R}, j = 1,\\ldots,\n","N]. (We implicitly assume that the $\\xi^i , i = 1,\\ldots, N$, are linearly independent; it follows that $W_N$ is an $N$ -dimensional subspace of $X^\\mathcal{N}$.) In the reduced-basis approach we look for an approximation $u_N (\\mu)$ to $u^\\mathcal{N} (\\mu)$ (which for our purposes here we presume is arbitrarily close to $u^e (\\mu))$ in $W_N$ ; in particular, we express $u_N (\\mu)$ as\n","\n$$\nu_N (\\mu) =\\sum_{i=1}^N u^j_N \\xi^j ;\n$$\n","we denote by latexmath:[%\n"," u_N (\\mu) \\in \\mathbb{R}^N] the coefficient vector $(u^1_N ,\\ldots, u^N_N )^T$. The premise — or hope — is that we should be able to accurately represent the solution at some new point in parameter space, $\\mu$, as an appropriate linear combination of solutions previously computed at a small number of points in parameter space (the $\\mu_i , i = 1,\\ldots, N ).$ But how do we find this appropriate linear combination? And how good is it? And how do we compute our approximation efficiently? The energy principle is crucial here (though more generally the weak form would suffice). To wit, we apply the classical Rayleigh-Ritz procedure to define\n","\n$$\n\\label{eq:3}\nu_N (\\mu) = \\mathop{\\mathrm{argmin}}_{wN \\in W_N} J(w_N );\n$$\n","alternatively we can apply Galerkin projection to obtain the equivalent statement\n","\n$$\n\\label{eq:4}\na(u_N (\\mu), v; \\mu) = \\ell(v),\\quad \\forall v \\in W_N .\n$$\n","The output can then be calculated from\n","\n$$\n\\label{eq:5}\n{T_{\\mathrm{root}}}_N (\\mu) = \\ell^O (u_N (\\mu)).\n$$\n","We now study this approximation in more detail.\n","a) Prove that, in the energy norm latexmath:[||| \\cdot ||| \\equiv (a(·, ·;\n","\\mu))^{1/2}] ,\n","\n$$\n|||u(\\mu) - u_N (\\mu)||| \\leq |||u(\\mu) - w_N |||, \\forall w_N \\in W_N .\n$$\n","This inequality indicates that out of all the possible choices of wN in the space $W_N$ , the reduced basis method defined above will choose the best one (in the energy norm). Equivalently, we can say that even if we knew the solution $u(\\mu)$, we would not be able to find a better approximation to $u(\\mu)$ in $W_N$ — in the energy norm — than latexmath:[u_N\n","(\\mu).]\n","b) Prove that\n","\n$$\nT_{\\mathrm{root}} (\\mu) - {T_{\\mathrm{root}}}_N (\\mu) = |||u(\\mu) - u_N (\\mu)|||^2.\n$$\n","c) Show that $u_N (\\mu)$ as defined in #eq:3-#eq:5 satisfies a set of $N \\times N$ linear equations,\n","\n$$\nA_N (\\mu) %\n u_N (\\mu) = %\n F_N ;\n$$\n","and that\n","\n$$\n{T_{\\mathrm{root}}}_N (\\mu) = %\n L^T_N %\n u_N (\\mu).\n$$\n","Give expressions for latexmath:[%\n"," A_N (\\mu) \\in \\mathbb{R}^{N\\times\n"," N}] in terms of latexmath:[%\n"," A^\\mathcal{N} (\\mu)] and latexmath:[Z,\n","%\n"," F^N \\in \\mathbb{R}^N] in terms of latexmath:[%\n"," F^\\mathcal{N}] and $Z,$ and latexmath:[%\n"," L^N \\in\n","\\mathbb{R}^N] in terms of latexmath:[%\n"," L^\\mathcal{N}] and $Z$; here $Z$ is an $\\mathcal{N} \\times N$ matrix, the $j$ th column of which is latexmath:[%\n"," u_N (\\mu^j)] (the nodal values of latexmath:[%\n"," u^\\mathcal{N}\n","(\\mu^j)]).\n","d) Show that the bilinear form $a(w, v; \\mu)$ can be decomposed as\n","\n$$\na(w, v; \\mu) = \\sum_{q=0}^Q \\theta^q (\\mu) a^q (w, v), \\forall w, v \\in X, \\forall \\mu \\in D,;\n$$\n","for $Q=6$ and give expressions for the $\\theta^q (\\mu)$ and the latexmath:[a^q\n","(w, v)]. Notice that the $aq (w, v)$ are not dependent on $\\mu$; the parameter dependence enters only through the functions latexmath:[\\theta^q\n","(\\mu), q = 1,\\ldots, Q.] Further show that\n","\n$$\n%\n A^\\mathcal{N} (\\mu) = \\sum_{q=1}^Q \\theta^q (\\mu) {%\n A^\\mathcal{N}}^q ,\n$$\n","and\n","\n$$\n%\n A^N (\\mu) = \\sum_{q=1}^Q \\theta^q (\\mu) %\n A^q_N ,\n$$\n","Give an expression for the latexmath:[{%\n"," A\\mathcal{N}}q] in terms of the nodal basis functions; and develop a formula for the latexmath:[%\n"," A^q_N] in terms of the latexmath:[{%\n"," A\\mathcal{N}}q] and $Z$.\n","e) The coercivity and continuity constants of the bilinear form for the continuous problem are denoted by $\\alpha^e (\\mu)$ and latexmath:[\\gamma^e\n","(\\mu)], respectively. We now assume that the basis function latexmath:[\\xi i ,\n","i = 1,\\ldots, N,] are orthonormalized. Show that the condition number of latexmath:[%\n"," A_N (\\mu)] is then bounded from above by latexmath:[\\gamma^e (\\mu)/\\alpha^e\n","(\\mu).]\n","f): Take into account the parameters $L$ and $t$ as parameters: the geometric transformation must be taken into account in the affine decomposition procedure.\n","*acknowledgment:* This problem set is based on a problem set in the class 16.920 Numerical Methods for Partial Differential Equations at MIT. We would like to thank Prof. A.T. Patera and Debbie Blanchard for providing the necessary material and the permission to use it.\n"],"metadata":{}}],"metadata":{"language_info":{"name":"python","version":"3.9.1"}},"nbformat":4,"nbformat_minor":4}
\ No newline at end of file
diff --git a/course-rom/env/antora.html b/course-rom/env/antora.html
index ed717fc..7f1b39c 100644
--- a/course-rom/env/antora.html
+++ b/course-rom/env/antora.html
@@ -145,9 +145,10 @@
Setting up the development Environment
diff --git a/course-rom/env/cmake.html b/course-rom/env/cmake.html
index 07a3c56..f82fbe2 100644
--- a/course-rom/env/cmake.html
+++ b/course-rom/env/cmake.html
@@ -145,9 +145,10 @@
Setting up the development Environment
diff --git a/course-rom/env/githubactions.html b/course-rom/env/githubactions.html
index 7f31dd2..3b99195 100644
--- a/course-rom/env/githubactions.html
+++ b/course-rom/env/githubactions.html
@@ -145,9 +145,10 @@
Setting up the development Environment
diff --git a/course-rom/env/jupyter.html b/course-rom/env/jupyter.html
index d39a008..2876484 100644
--- a/course-rom/env/jupyter.html
+++ b/course-rom/env/jupyter.html
@@ -145,9 +145,10 @@
Setting up the development Environment
diff --git a/course-rom/env/rename.html b/course-rom/env/rename.html
index b4d3992..b3e25e3 100644
--- a/course-rom/env/rename.html
+++ b/course-rom/env/rename.html
@@ -145,9 +145,10 @@
Setting up the development Environment
diff --git a/course-rom/env/vscode.html b/course-rom/env/vscode.html
index 9af16d7..c9c4abd 100644
--- a/course-rom/env/vscode.html
+++ b/course-rom/env/vscode.html
@@ -145,9 +145,10 @@
Setting up the development Environment
diff --git a/course-rom/homework/problem-set-1.html b/course-rom/homework-2023/problem-set-1.html
similarity index 98%
rename from course-rom/homework/problem-set-1.html
rename to course-rom/homework-2023/problem-set-1.html
index a7f2502..5e5647a 100644
--- a/course-rom/homework/problem-set-1.html
+++ b/course-rom/homework-2023/problem-set-1.html
@@ -4,7 +4,7 @@
Problem Set 1: RB for Linear Affine Elliptic Problems :: Course Rom
-
+
@@ -145,9 +145,10 @@
-
.ipynb
diff --git a/course-rom/homework/problem-set-2.html b/course-rom/homework-2023/problem-set-2.html
similarity index 98%
rename from course-rom/homework/problem-set-2.html
rename to course-rom/homework-2023/problem-set-2.html
index a360f13..e0ceb59 100644
--- a/course-rom/homework/problem-set-2.html
+++ b/course-rom/homework-2023/problem-set-2.html
@@ -4,7 +4,7 @@
Problem Set 2: RB for Linear Affine Elliptic Problems :: Course Rom
-
+
@@ -147,9 +147,10 @@
Problem Set 1: RB for Linear Affine Elliptic Problems
+
+
+
+
We consider the problem of designing a thermal fin to effectively remove heat from a surface. The two-dimensional fin, shown in Figure 1, consists of a vertical central “post” and four horizontal “subfins”; the fin conducts heat from a prescribed uniform flux “source” at the root, \(\Gamma_{\mathrm{root}}\) , through the large-surface-area subfins to surrounding flowing air. The fin is characterized by a five-component parameter vector, or “input,” \(\mu_
+= (\mu_1 , \mu_2, \ldots, \mu_5 )\),where \(\mu_i = k^i , i = 1, \ldots
+, 4\), and \(\mu_5 = \mathrm{Bi}\); \(\mu\) may take on any value in a specified design set \(D \subset \mathbb{R}^5\).
+
+
+
+
+
+
Figure 1. Thermal fin
+
+
+
Here \(k^i\) is the thermal conductivity of the ith subfin (normalized relative to the post conductivity \(k^0 \equiv 1\)); and \(\mathrm{Bi}\) is the Biot number, a nondimensional heat transfer coefficient reflecting convective transport to the air at the fin surfaces (larger \(\mathrm{Bi}\) means better heat transfer). For example, suppose we choose a thermal fin with \(k^1 = 0.4, k^2 = 0.6, k^3 = 0.8, k^4 = 1.2\), and \(\mathrm{Bi} = 0.1\); for this particular configuration \(\mu = \{0.4, 0.6, 0.8, 1.2, 0.1\}\), which corresponds to a single point in the set of all possible configurations D (the parameter or design set). The post is of width unity and height four; the subfins are of fixed thickness \(t = 0.25\) and length \(L = 2.5\).
+
+
+
We are interested in the design of this thermal fin, and we thus need to look at certain outputs or cost-functionals of the temperature as a function of \(\mu\). We choose for our output \(T_{\mathrm{root}}\), the average steady-state temperature of the fin root normalized by the prescribed heat flux into the fin root. The particular output chosen relates directly to the cooling efficiency of the fin — lower values of \(T_{\mathrm{root}}\) imply better thermal performance. The steady–state temperature distribution within the fin, \(u(\mu)\), is governed by the elliptic partial differential equation
+
+
+
+\[\label{eq:1}
+-k^i \Delta u^i = 0 \text{ in } \Omega^i , i = 0, \ldots, 4,\]
+
+
+
+
where \(\Delta\) is the Laplacian operator, and \(u_i\) refers to the restriction of \(u \text{ to } \Omega^i\) . Here \(\Omega^i\) is the region of the fin with conductivity \(k^i , i = 0,\ldots, 4: \Omega^0\) is thus the central post, and \(\Omega^i , i = 1,\ldots, 4\), corresponds to the four subfins. The entire fin domain is denoted \(\Omega (\bar{\Omega} = \cup_{i=0}^4 \bar{\Omega}^i )\); the boundary \(\Omega\) is denoted \(\Gamma\). We must also ensure continuity of temperature and heat flux at the conductivity– discontinuity interfaces \(\Gamma^i_{int} \equiv \partial\Omega^0 \cap \partial\Omega^i , i = 1,\ldots, 4\), where \(\partial\Omega^i\) denotes the boundary of \(\Omega^i\), we have on \(\Gamma^i_{int} i = 1,\ldots, 4\) :
which models the heat source; and a Robin boundary condition
+
+
+
+\[-k^i (\nabla u^i \cdot n^i ) = \mathrm{Bi} u^i \text{ on } \Gamma^i_{ext} , i = 0,\ldots, 4,\]
+
+
+
+
which models the convective heat losses. Here \(\Gamma^i_{ext}\) is that part of the boundary of \(\Omega^i\) exposed to the flowing fluid; note that \(\cup_{i=0}^4 \Gamma^i_{ext} = \Gamma\backslash\Gamma_{\mathrm{root}}\). The average temperature at the root, \(T_{\mathrm{root}} (\mu)\), can then be expressed as \(\ell^O(u(\mu))\), where
(recall \(\Gamma_{\mathrm{root}}\) is of length unity). Note that \(\ell(v) = \ell^O(v)\) for this problem.
+
+
+
+
+
1. Part 1 - Finite Element Approximation
+
+
+
We saw in class that the reduced basis approximation is based on a “truth” finite element approximation of the exact (or analytic) problem statement. To begin, we have to show that the exact problem described above does indeed satisfy the affine parameter dependence and thus fits into the framework shown in class.
+
+
+
a) Show that \(u^e (\mu) \in X^e \equiv H^1 (\Omega)\) satisfies the weak form
+\[\begin{aligned}
+a(w, v; \mu) &=\sum_{i=0}^4 k^i \int_{\Omega^i} \nabla w \cdot \nabla v dA
++ \mathrm{Bi}\int_{\Gamma\backslash\Gamma_{\mathrm{root}}} w v dS,\\
+\ell(v) &= \int_{\Gamma_{\mathrm{root}}} v
+\end{aligned}\]
+
+
+
+
b) Show that \(u^e (\mu)\) is the argument that minimizes
+
+
+
+\[\label{eq:2}
+J(w) = \frac{1}{2}\sum_{i=0}^4 k^i \int_{\Omega_i} \nabla w \cdot
+\nabla w dA + \frac{\mathrm{Bi}}{2}
+\int_{\Gamma\backslash\Gamma_{\mathrm{root}}} w^2 dS -
+\int_{\Gamma_{\mathrm{root}}} w dS\]
+
where \(%
+ A^{\mathcal{N}} \in
+\mathbb{R}^{\mathcal{N}\times\mathcal{N}} ,
+%
+ u^{\mathcal{N}} \in \mathbb{R}^{\mathcal{N}} ,
+%
+ F^{\mathcal{N}} \in \mathbb{R}^{\mathcal{N}} , and
+%
+ L^{\mathcal{N}} \in \mathbb{R}^{\mathcal{N}}\); here \(\mathcal{N}\) is the dimension of the finite element space \(X^{\mathcal{N}}\), which (given our natural boundary conditions) is equal to the number of nodes in \(\mathcal{T}_h\).
+
+
+
+
+
2. Part 2 - Reduced-Basis Approximation
+
+
+
In general, the dimension of the finite element space, \(\mathop{\mathrm{dim}}X =
+\mathcal{N}\), will be quite large (in particular if we were to treat the more realistic three-dimensional fin problem), and thus the solution of \(%
+ A^{\mathcal{N}} %
+ u^{\mathcal{N}} (\mu) = %
+ F^\mathcal{N}\) can be quite expensive. We thus investigate the reduced-basis methods that allow us to accurately and very rapidly predict \(T_{\mathrm{root}} (\mu)\) in the limit of many evaluations — that is, at many different values of \(\mu\) — which is precisely the “limit of interest” in design and optimization studies. To derive the reduced-basis approximation we shall exploit the energy principle,
To begin, we introduce a sample in parameter space, \(S_N = \{\mu_1 ,
+ \mu_2 ,\ldots, \mu_N \}\) with \(N \ll \mathcal{N}\). Each \(\mu_i , i =
+1,\ldots, N\) , belongs in the parameter set \(\mathcal{D}\). For our parameter set we choose \(\mathcal{D} = [0.1, 10.0\)^4 \times [0.01, 1.0]], that is, \(0.1 \leq k^i \leq 10.0, i = 1,\ldots, 4\) for the conductivities, and \(0.01
+\leq \mathrm{Bi} \leq 1.0\) for the Biot number. We then introduce the reduced-basis space as
where \(u_N (\mu_i )\) is the finite-element solution for \(\mu =
+\mu_i\).
+
+
+
To simplify the notation, we define \(\xi^i \in X\) as \(\xi^i = u_N
+(\mu_i ), i = 1,\ldots, N\); we can then write \(W_N = span\{\xi^i , i =
+ 1,\ldots, N \}\). Recall that \(W_N = \mathop{\mathrm{span}}\{\xi^i , i = 1,\ldots, N \}\) means that \(W_N\) consists of all functions in \(X\) that can be expressed as a linear combination of the \(\xi^i\) ; that is, any member \(v_N\) of \(W_N\) can be represented as
+
+
+
+\[v_N =\sum_{i=1}^N \beta^j \xi^j ,\]
+
+
+
+
for some unique choice of \(\beta^j \in \mathbb{R}, j = 1,\ldots,
+N\). (We implicitly assume that the \(\xi^i , i = 1,\ldots, N\), are linearly independent; it follows that \(W_N\) is an \(N\) -dimensional subspace of \(X^\mathcal{N}\).) In the reduced-basis approach we look for an approximation \(u_N (\mu)\) to \(u^\mathcal{N} (\mu)\) (which for our purposes here we presume is arbitrarily close to \(u^e (\mu))\) in \(W_N\) ; in particular, we express \(u_N (\mu)\) as
+
+
+
+\[u_N (\mu) =\sum_{i=1}^N u^j_N \xi^j ;\]
+
+
+
+
we denote by \(%
+ u_N (\mu) \in \mathbb{R}^N\) the coefficient vector \((u^1_N ,\ldots, u^N_N )^T\). The premise — or hope — is that we should be able to accurately represent the solution at some new point in parameter space, \(\mu\), as an appropriate linear combination of solutions previously computed at a small number of points in parameter space (the \(\mu_i , i = 1,\ldots, N ).\) But how do we find this appropriate linear combination? And how good is it? And how do we compute our approximation efficiently? The energy principle is crucial here (though more generally the weak form would suffice). To wit, we apply the classical Rayleigh-Ritz procedure to define
This inequality indicates that out of all the possible choices of wN in the space \(W_N\) , the reduced basis method defined above will choose the “best one” (in the energy norm). Equivalently, we can say that even if we knew the solution \(u(\mu)\), we would not be able to find a better approximation to \(u(\mu)\) in \(W_N\) — in the energy norm — than \(u_N
+(\mu).\)
Give expressions for \(%
+ A_N (\mu) \in \mathbb{R}^{N\times
+ N}\) in terms of \(%
+ A^\mathcal{N} (\mu)\) and \(Z,
+%
+ F^N \in \mathbb{R}^N\) in terms of \(%
+ F^\mathcal{N}\) and \(Z,\) and \(%
+ L^N \in
+\mathbb{R}^N\) in terms of \(%
+ L^\mathcal{N}\) and \(Z\); here \(Z\) is an \(\mathcal{N} \times N\) matrix, the \(j\) th column of which is \(%
+ u_N (\mu^j)\) (the nodal values of \(%
+ u^\mathcal{N}
+(\mu^j)\)).
+
+
+
d) Show that the bilinear form \(a(w, v; \mu)\) can be decomposed as
for \(Q=6\) and give expressions for the \(\theta^q (\mu)\) and the \(a^q
+(w, v)\). Notice that the \(aq (w, v)\) are not dependent on \(\mu\); the parameter dependence enters only through the functions \(\theta^q
+(\mu), q = 1,\ldots, Q.\) Further show that
Give an expression for the \({%
+ A^\mathcal{N}}^q\) in terms of the nodal basis functions; and develop a formula for the \(%
+ A^q_N\) in terms of the \({%
+ A^\mathcal{N}}^q\) and \(Z\).
+
+
+
e) The coercivity and continuity constants of the bilinear form for the continuous problem are denoted by \(\alpha^e (\mu)\) and \(\gamma^e
+(\mu)\), respectively. We now assume that the basis function \(\xi i ,
+i = 1,\ldots, N,\) are orthonormalized. Show that the condition number of \(%
+ A_N (\mu)\) is then bounded from above by \(\gamma^e (\mu)/\alpha^e
+(\mu).\)
+
+
+
f): Take into account the parameters \(L\) and \(t\) as parameters: the geometric transformation must be taken into account in the affine decomposition procedure.
+
+
+
acknowledgment: This problem set is based on a problem set in the class “16.920 Numerical Methods for Partial Differential Equations” at MIT. We would like to thank Prof. A.T. Patera and Debbie Blanchard for providing the necessary material and the permission to use it.
Problem Set 2: RB for Linear Affine Elliptic Problems
+
+
1. Problem Statement — Design of a Thermal Fin
+
+
+
We consider the problem of designing a thermal fin described in Problem Set 1. In PS1 we looked at some thoeretical issues (weak formulation and optimization formulation, convergence of the reduced basis approximation) and derived the necessary reduced basis quantities, i.e., expressions for \(A_N ( \mu )\), \(F_N\) , and \(L_N\) . This problem set is devoted to implementing the reduced basis approximation and solving a simple design problem.
+
+
+
+
+
2. Reduced Basis Approximation
+
+
+
The point of departure for the reduced basis approximation is a high – dimensional finite element “truth” discretization. In the offline stage we require the finite element solution to build the reduced basis and we thus also need the FE matrices. In this problem set we skip the FE assembly step and provide all of the necessary data for use in Python (see Appendix 1).
+
+
+
We saw in class that the reduced basis solution \(u_N ( \mu ) \in \mathbb{R}^N\) satisfies the set of \(N\times N\) linear equations,
+
+
+
+\[ A_N( \mu )u_N( \mu ) = F_N;\]
+
+
+
+
and that the outputis given by
+
+
+
+\[ {T_{root}}_N ( \mu ) = L^T_N u_N ( \mu ).\]
+
+
+
+
We derived expressions for \(A_N( \mu ) \in \mathbb{R}^{N\times N}\) in terms of \(A_N( \mu )\) and \(Z\), \(F_N \in \mathbb{R}^N\) in terms of \(F_N\) and \(Z\), and \(L_N \in \mathbb{R}^N\) in terms of \(L_N\) and \(Z\); here \(Z\) is an \(\mathcal{N} \times N\) matrix, the jth column of which is \(u_N ( \mu j )\) (the nodal values of \(u_N ( \mu j ))\). Finally, it follows from affine parameter dependence that \(A_N ( \mu )\) can be expressed as
Evaluate the output \({T_{root}}_N ( \mu )\) from 1.2).
+
+
+
+
+
+
+
+
1 The idea is that the offline stage is done only once, generating a small datafile with the \(A^q_N , q = 1,\ldots,Q\), \(F_N\), and \(L_N\); the on-line stage then accesses this datafile to provide real-time response to new \(\mu\) queries. For the required off-line finite element calculations in this and the following questions, you should first use the coarse triangulation \(\mathcal{T}_{h,\mathrm{coarse}}\).
+
+
+
+
+
Show that the operation count for the on-line stage of your code is independent of \(\mathcal{N}\) . In particular show that the operation count (number of floating-point operations) for the on-line stage, for each new \(\mu\) of interest, can be expressed as
for \(c_1, c_2, c_3, \gamma_1, \gamma_2,\) and \(\gamma_3\) independent of \(N\). Give values for the constants \(c_1, c_2, c_3, \gamma_1, \gamma_2,\) and \(\gamma_3\).
+
+
+
+
+
We first consider a one parameter (\(P = 1\)) problem. To this end, we keep the Biot number fixed at \(Bi = 0.1\) and assume that the conductivities of all fins are equivalent, i.e., \(k_1 = k_2 = k_3 = k_4\), but are allowed to vary between \(0.1\) and \(10\) – we thus have \(\mu \in D =
+[0.1, 10].\) The sample set \(S_N\) for \(N_{max} = 8\) is given the log equidistributed sampling.
+
+
+
Generate the reduced basis “matrix” \(Z\) and all necessary reduced basis quantities. You have two options: you can use the solution "snapshots" directly in \(Z\) or perform a Gram-Schmidt orthonormalization to construct \(Z\) (Note that you require the \(X\) – inner product to perform Gram-Schmidt; here, we use \((\cdot, \cdot)_X = a(\cdot, \cdot; \mu )\), where \(\mu = 1\) – all conductivities are \(1\) and the Biot number is \(0.1\)). Calculate the condition number of \(A_N ( \mu )\) for \(N = 8\) and for \(\mu = 1\) and \(\mu = 10\) with and without Gram – Schmidt orthonormalization. What do you observe? Solve the reduced basis approximation (where you use the snapshots directly in \(Z\)) for \(\mu_1 = 0.1\) and \(N = 8\). What is \(u_N( \mu_1)\)? How do you expect \(u_N( \mu_2)\) to look like for \(\mu_2
+= 10.0\)? What about \(\mu_3 = 1.0975\)? Solve the Gram – Schmidt orthonormalized reduced basis approximation for \(\mu_1 = 0.1\) and \(\mu
+2 = 10\) for \(N = 8\). What do you observe? Can you justify the result? For the remaining questions you should use the Gram – Schmidt orthonormalized reduced basis approximation.
+
+
+
+
Verify that, for \(\mu = 1.5\) (recall that Biot is still fixed at \(0.1\)) and \(N = 8\), the value of the output is \({T_{root}}_N ( \mu ) = 1.61\) up to 2 digits.
+
+
+
We next introduce a regular test sample, \(\Xi_{test} \subset D\), of size \(ntest = 100\) (in Python you can simply use linspace(0.1, 10, 100) to generate \(\Xi_{test}\)). Plot the convergence of the maximum relative error in the energy norm \(\max_{\mu \in\Xi_{test}} |||u( \mu ) -
+u_N ( \mu )|||_\mu /|||u( \mu )|||_\mu\) and the maximum relative output error max \(\mu \in\Xi_{test} |{T_{root}}( \mu ) - {T_{root}} N( \mu
+)|/{T_{root}}( \mu )\) as a function of \(N\) (use the Python command semilogy for plotting).
+
+
+
Compare the average CPU time over the test sample required to solve the reduced basis online stage with direct solution of the FE approximation as a function of \(N\).
+
+
+
What value of \(N\) do you require to achieve a relative accuracy in the output of 1%. What savings in terms of CPU time does this % correspond to?
+
+
+
Solve problems b) 3. to 5. using the medium and fine FE triangulation. Is the dependence on \(\mathcal{N}\) as you would anticipate?
+
+
+
+
+
+
We now consider another one parameter \((P = 1)\) problem. This time, we assume that the conductivities are fixed at \(\{k_1,k_2,k_3,k_4\} = \{0.4,0.6,0.8,1.2\}\), and that only the Biot number, \(Bi\), is allowed to vary from \(0.01\) to \(1\). The sample set \(S_N\) for \(N_{max} = 11\) is given by log equidistributed sampling. Generate an orthonormal \(Z\) from the sample set using the medium triangulation.
+
+
+
+
Verify that, for \(\mu_0 = {0.4, 0.6, 0.8, 1.2, 0.15}\), i.e. \(Bi = 0.15\), the value of the output is \({T_{root}}_N ( \mu 0) =1.53\).
+
+
+
We next introduce a regular test sample, \(\Xi_{test} \subset D\), of size \(ntest =100\) (in Python you can simply use linspace(0.01, 1, 100) to generate \(\Xi_{test}\)). Plot the convergence of the maximum relative error in the energy norm \(\max_{\mu \in\Xi_{test}} |||u( \mu ) - u_N ( \mu )|||_\mu /|||u( \mu
+)|||_\mu\) and the maximum relative output error \(\max_{\mu \in\Xi_{test}} |{T_{root}}( \mu ) - {T_{root}}_N( \mu )|/{T_{root}}(
+ \mu )\) as a function of \(N\) (use the Python command semilogy for plotting).
+
+
+
The Biot number is directly related to the cooling method; higher cooling rates (higher \(Bi\)) imply lower (better) \({T_{root}}\) but also higher (worse) initial and operational costs. We can thus define (say) a total cost function as
+
+
+\[C(Bi) = Bi + {T_{root}}(Bi),\]
+
+
+
+
minimization of which yields an optimal solution. Apply your (online) reduced – basis approx – imation for \({T_{root}}_N\) (that is, replace \({T_{root}}(Bi)\) in (above) with \({T_{root}}_N (Bi))\) to find the optimal \(Bi.\) Any (simple) optimization procedure suffices for the minimization.
+
+
+
+
+
+
+
We consider now a two parameter \((P = 2)\) problem where the conductivities are assumed to be equivalent, i.e., \(k_1 = k_2 = k_3 = k_4\), but are allowed to vary between \(0.1\) and \(10\); and the Biot number, \(Bi\), is allowed to vary from \(0.01\) to \(1\). The sample set \(S_N\) for \(N_{max} = 46\) is given by the log random sampling. Generate an orthonormal \(Z\) from the sample set using the coarse triangulation.
+
+
+
We next introduce a regular grid, \(\Xi_{test} \subset D\), of size \(ntest = 400\) (a regular \(20 \times 20\) grid). Plot the convergence of the maximum relative error in the energy norm \(\max_{\mu \in\Xi_{test}} |||u( \mu ) - u_N ( \mu )|||_\mu /|||u( \mu
+)|||_\mu\) and the maximum relative output error \(max_{\mu \in \Xi_{test}} |{T_{root}}( \mu ) - {T_{root}}_N( \mu )|/{T_{root}}( \mu)\) as a function of \(N\).
+
+
+
We now consider the POD method and we wish to compare it with the Greedy approximation. To this end, we sample log randomly the parameter space (\(P=2\)) and take \(n_{\mathrm{train}}=100\) samples. Build the POD approximation using these samples as training set and compare the results with the Greedy approximation. Compute the RIC and the dimension of the POD space (\(N\)) such that the RIC is \(99\%\) of the total energy. Plot the POD and Greedy convergence of the maximum relative error in the energy norm \(\max_{\mu \in\Xi_{test}} |||u( \mu ) - u_N ( \mu )|||_\mu /|||u( \mu
+)|||_\mu\) and the maximum relative output error \(max_{\mu \in \Xi_{test}} |{T_{root}}( \mu ) - {T_{root}}_N( \mu )|/{T_{root}}( \mu
+)\) as a function of \(N\).
+
+
+
Implement the parametrisation with respect to \(L\) and \(t\). The reference geometry is the one given by the .geo file and the corresponding \(\hat{L}\) and \(\hat{t}\). Plot the mean temperature \({T_{root}}( \mu )\) as a function \(t \in [0.1,0.5]\) and the other parameters set to \(k_i=0.1, L=2.5, Bi=0.1\).
+
+
+
+
+
+
+
3. A Posteriori Error Bounds, Greedy Sampling Procedure
+
+
+
We consider again the problem of designing a thermal fin of Problem Set 1 and 2. Given the reduced basis approximation implemented in PS2, we turn to implementing the associated a posteriori error estimation procedures developed in the lecture. The second half of this problem set is devoted to implementing the greedy sampling procedure. We will consider the following two cases:
+
+
+
+
Case I (\(P=1\))
+
+
We keep the Biot number fixed at \(Bi = 0.1\) and assume that the conductivities of all fins are equivalent, i.e., \(k = k_1 = k_2 = k_3 = k_4\), but are allowed to vary between 0.1 and 10 — we thus have \(\mu \in D = [0.1, 10\).] For this \(P = 1\) case we define the \(X\)-inner product \((\cdot, \cdot)_X = a(\cdot, \cdot; \bar{\mu}),\) where \(\bar{\mu} = 1.\)
+
+
+
+
+
We also define the parameter grids \(G^{\mathrm{lin}}_{[ \mu_{min} , \mu_{max} ;10]}\) and \(G^{\mathrm{ln}}_{[ \mu_{min} , \mu_{max} ;10]}\). The former grid is equi-spaced in \(\mu\), the latter grid is equi-spaced in \(ln(\mu)\) — often advantageous within the reduced basis context. More generally, the “log” spacing represents equal relative increments, and thus represents better coverage for parameters that vary over a large range. For the \(P = 2\) case we can then define tensor grids over \(\mathcal{D}\), \(\Xi^{\mathrm{log}}_M \subset D \subset \mathbb{R}^2\) , as
note \(\Xi^{log}_M\) contains \(M^2\) points; a similar definition applies to \(\Xi^{lin}_M\). We also define a particular test grid (biased neither towards “log” nor “lin”)
Given the coercivity lower bound, we can now turn to implementing the a posteriori error bounds. Note that, in principle, there is an online-inefficient and an online-efficient way to evaluate these error bounds. We first consider the latter: From the lecture we know that the energy norm a posteriori error bound is given by
+\[r(v; \mu) = f (v; \mu) - a(u_N (\mu), v; \mu),\quad \forall v \in X.\]
+
+
+
+
For any new \(\mu\) and associated reduced basis solution, \(u_N (\mu),\) we can now directly calculate \(\hat{e}(\mu)\) from 2.2 and 2.3, evaluate the norm \(\|\hat{e}(\mu)||_X\) and — given \(\alpha_{LB} (\mu)\) — obtain \(\Delta^{en}_N (\mu)\) from 2.1. Although this approach is online-inefficient because the computational cost depends on \(O(\mathcal{N})\), it is nevertheless useful as a means to test your offline-online computational decomposition. We will consider Case I and Case II in the sequel. Note that you should only require one code to handle both cases, i.e., Case I is a specialization of Case II by keeping one of the parameters fixed. Also, when using you should try to replace loops by matrix-vector products as much as possible (e.g. try to write the nested loop over \(N\) when summing up the contributions of the \(\|\hat{e}(\mu)\|_X\) norm as a vector-matrix-vector product — the nested loop over \(Q_a\) is less critical).
+
+
+
We first consider Case I.
+To answer this question you should use the sample set \(S_N\) provided for PS2 (RB_sample.sample1), orthonormalize the basis functions, and use the medium grid.
+
+
+
+
+
Implement an offline/online version of the a posteriori error bound calculation (not using the affine decomposition) shown in the lecture (this is inefficient). Compute the direct calculation for the error bound, \(\Delta^{en}_N (\mu)\), for all \(N (1 \leq N \leq 8)\) and (say) \(5\) parameter values randomly distributed in \(\mathcal{D}.\)
+
+
+
Calculate \(\eta^{en}_{\min,N},\eta^{en}_{\max,N}\) and \(\eta^{en}_{ave,N}\) the minimum, maximum, and average effectivity \(\eta^{en}_N(\mu)\) over \(\Xi test = G^{lin}[ \mu_{min} , \mu_{max} ;50\) \cup G^{ln}[ \mu_{min} , \mu_{max} ;50\]], respectively (note that \(\Xi^{test}\) is of size 20 since \(P = 1\)).
+
+
+
+
+
Present the results in a table for all \(N\) . Is the minimum effectivity greater than unity? How does the maximum effectivity compare with your theoretical upper bound for the effectivity? (Note you should exclude from the min/max/mean-operation all points in \(\Xi^{test}\) for which \(\|u(\mu) - u_N (\mu)\|_X\) is less than (say) \(10e-11\) .)
+
+
+
+
+
4. Parablic problem
+
+
+
4.1. Thermal Fin Problem
+
+
Our problem of interest is the thermal fin discussed in the previous problem sets, but now we consider the time-dependent case. We assume that the thermal fin is initially at zero (non-dimensionalized) temperature and a heat flux is then applied to the root. The output of interest is the average temperature of the fin. We directly consider the truth approximation. To this end, we divide the time interval, \(I = (0,t_f\)], into K subintervals of equal length \(\Delta = \frac{t_f}{K}\), and define \(t_k = k\Delta t, 0 \leq k \leq K\). We shall consider Euler-Backward for the time integration. We also recall the truth finite element approximation space \(X \subset X^e\).
+
+
+
Our truth problem statement is then: given a parameter \(\mu \in D\), we evaluate the output:
+
+
+
\(s^k(\mu) = l(u^k(\mu)), 1 \leq k \leq K\)
+
+
+
where the field variable \(u^k(\mu) \in X, 1 \leq k \leq K\), satisfies:
+
+
+
\(m\left( \frac{u^k(\mu)-u^{k-1}(\mu)}{\Delta t}, v \right) + a( u^k(\mu), v; \mu ) = f(v)g(t^k), \forall v \in X\)
+
+
+
with initial condition \(u(t_0; \mu) = u_0 = 0\).
+
+
+
Here:
+- The bilinear form \(a\) is defined as in Problem Set 1.
+- The linear form \(f\) is given by \(f(v) = \int_{\Gamma_{root}}v\).
+- The linear form \(l\) is given by \(l(v) = \int_\Omega v\).
+- The bilinear form \(m\) is given by:
+
+
+
\(m(u,v) = \int_\Omega uv, \forall u, v \in X\)
+
+
+
and \(g(t_k)\) denotes the "control input" at time \(t = t_k\). Note that \(m\) and \(l\), \(f\) are parameter-independent.
+
+
+
We consider the following special case: We assume that the conductivities of all fins are equivalent and fixed at \(k_i = 1, i = 1,...,4,\) and that the Biot number is allowed to vary between 0.01 and 1. We thus have \(\mu \equiv Bi \in D = [0.01, 1\)]. We consider the time interval \(I = (0, 10\)] with a discrete timestep \(\Delta t = 0.1\) and thus \(K = 100\).
+
+
+
To begin, you should download and unpack the zip file PS4_matlab.zip. You will find the file FE_matrix_mass.mat which contains a struct, FE_matrix_mass, with the mass matrices for the fine, medium, and coarse triangulations used before. To generate the output vector \(L\) you can simply postmultiply the corresponding mass matrix with a vector containing all 1s. From the previous problem sets you already have the required finite element forcing vector F and the finite element stiffness matrix \(A\) (and the \(A_q\)). In the sequel, you should use the medium triangulation.
+
+
+
+
4.2. Part 1 - Reduced Basis Approximation
+
+
We first generate a reduced basis approximation by choosing a basis from scratch. To this end, we use \(g(t_k) = \delta_{1k}, 1 \leq k \leq 100\) (unit impulse input) and set:
i.e., our reduced basis space \(X_N\) is spanned by the solution \(u^k(\mu)\) at several parameter-time pairs. We then orthonormalize \(X_N\) using Gram-Schmidt.
+
+
+
4.2.1. Tasks
+
+
+
+
Write an offline-online code in MATLAB for the reduced basis approximation (use LU decomposition for the truth and reduced basis time integration).
+
+
+
+
+
+
# generate RB code
+
+
+
+
+
+
Plot the outputs \(s^k(\mu)\), \(s^k_N(\mu)\), and the error \(s^k(\mu) - s^k_N(\mu)\) as a function of time for \(g(t_k) = 1 - \cos(t_k)\) and \(\mu = 0.05\).
+
+
+
Plot \(|||u^k(\mu)|||\), \(|||u^k_N(\mu)|||\), and the error \(|||u^k(\mu) - u^k_N(\mu)|||\) as a function of time for \(g(t_k) = 1 - \cos(t_k)\) and \(\mu = 0.05\).
+
+
+
+
+
+
+
4.3. Part 2 - A Posterior Error Estimation
+
+
The problem statement fits in the framework introduced in the lecture.
+
+
+
+
4.4. Tasks
+
+
+
+
Similar to the elliptic case, derive and implement an offline-online version for the direct calculation (no offline/online calculation) of the energy norm a posteriori error bound for the primal variable by extending your code from the elliptic case.
+
+
+
+
+
+
# python code
+
+
+
+
+
+
Compute the direct calculation of the error bound for 10 random parameter values in D. You can perform this over time (better) or compare the values at the final time.
+
+
+
+
+
+
4.5. Part 3 - Sampling Procedure
+
+
Our reduced basis space from Part 1 is less than optimal. Given your offline-online decomposition for the reduced basis approximation from Part 1 and associated a posteriori error estimation from Part 2, we can now pick a much more optimal basis.
+
+
+
4.5.1. Tasks
+
+
+
+
Apply the POD-Greedy algorithm with \(\Xi_{train} = G^{ln}_{[0.01,1;100\)}, \varepsilon^{tol,min}=1e-6,] and \(\mu^{*}_0 = 0.01\). Here, we also use the impulse input \(g(t_k) = \delta_{1k}, 1 \leq k \leq 100\).
+
+
+
+
Determine \(N_{max}\) to achieve the desired accuracy.
+
+
+
Plot \(\Delta^{max}_N = \Delta^K_N(\mu^{*})/|||u^k(\mu^*)|||\) as a function of \(N\).
+
+
+
Plot the outputs \(s^k(\mu)\), \(s^k_N(\mu)\), and the simple error bound \(\Delta^s_N(t_k; \mu)\).
+
+
+
+
+
+
+
+
+
+
+
+
5. Appendix
+
+
+
Here’s the converted content in AsciiDoc format using stem for LaTeX math:
+
+
+
+
+
6. Implementation of the Reduced Basis Method
+
+
+
For the implementation of the reduced basis method, the finite element matrices for three possible triangulations of the fin problem are provided. To obtain the required Python data, download the file PS2_Python.zip from the course website and unzip it. There are three .mat files: FE_matrix.mat contains the FE matrices, FE_grid.mat contains the triangulation data, and RB_sample.mat contains the samples you should use initially (later on, you will generate samples yourselves using a greedy procedure). To load the FE matrices in the Python workspace:
+
+
+
+
load FE_matrix
+
+
+
+
This creates one variable named FE_matrix with three fields: coarse, medium, and fine. Each of these fields contains a cell array Ahq of size \(6 \times 1\) and the load vector Fh. Each cell of Ahq contains the parameter-independent FE matrix \(A_N^q\), where \(q = 1, ..., 6\). Here:
+
+
+
+
+
\(q = 1, ..., 4\): Corresponds to the "submatrices" of fins \(1, ..., 4\) with conductivities \(k_i, i = 1, ..., 4\), respectively.
+
+
+
\(q = 5\): Corresponds to the "submatrix" of the central post with conductivity \(k_0 = 1\).
+
+
+
\(q = 6\): Corresponds to the "submatrix" of the line integral over the "surface" of the fin (without \(\Gamma_{root}\)).
+
+
+
+
+
To load the reduced basis samples \(S_N\) in the Python workspace:
+
+
+
+
load RB_sample
+
+
+
+
This creates one variable named RB_sample with fields sample1, sample2, and sample3, corresponding to the two \(P = 1\) cases and the \(P = 2\) case described in the problem statement.
+
+
+
Note that you require the triangulation only for plotting the FE solution (see below). The following detailed information about the triangulation is just included to give you an impression concerning the data required if you would like to set up the FE matrices from scratch. To load the triangulation data in the Python workspace:
+
+
+
+
load FE_grid
+
+
+
+
This creates one variable FE_grid with three fields: coarse, medium, and fine. Each of these fields is a different triangulation \(\mathcal{T}_h\) for the fin problem. More specifically:
+
+
+
+
+
coarse defines \(\mathcal{T}_{h,coarse}\) with 1333 nodes and 2095 elements.
+
+
+
medium defines \(\mathcal{T}_{h,medium}\) with 4760 nodes and 8380 elements.
+
+
+
fine defines \(\mathcal{T}_{h,fine}\) with 17889 nodes and 33520 elements.
+
+
+
+
+
Each of these variables is of type struct, with four different fields:
Description of the fields (assume that we are using the coarse triangulation):
+
+
+
+
+
nodes: The number of nodes in the triangulation.
+
+
+
coor: Two-dimensional matrix with size (nodes × 2), where each row \(i\) has the \(x\) and \(y\) coordinates for node \(i\). For example, the location of node 49 can be determined by two coordinates. The coordinate in the \(x\)-direction would be coarse.coor(49,1) and in the \(y\)-direction coarse.coor(49,2).
+
+
+
elements: The number of elements in the triangulation.
+
+
+
theta: The adjacency matrix \(\Theta(k, \alpha)\) which defines the local-to-global mapping required in the elemental assembly procedure. Since we have regions with different physical properties, for each region a separate adjacency matrix is provided. The regions considered are:
+
+
+
Region 1: Domain \(\Omega_1, \Theta_1(k, \alpha) = coarse.theta{1}\)
+
+
+
Region 2: Domain \(\Omega_2, \Theta_2(k, \alpha) = coarse.theta{2}\)
+
+
+
Region 3: Domain \(\Omega_3, \Theta_3(k, \alpha) = coarse.theta{3}\)
+
+
+
Region 4: Domain \(\Omega_4, \Theta_4(k, \alpha) = coarse.theta{4}\)
+
+
+
Region 5: Domain \(\Omega_0, \Theta_5(k, \alpha) = coarse.theta{5}\)
+
+
+
+
+
For each of these regions \(i\), the index \(k\) varies in the range \(k \in \{1, ..., n_i\}\), where \(n_i\) is the number of elements in region \(i\). For example, element 12 in region 3 has the global nodes \(\nu_1 = coarse.theta{3}(12,1)\), \(\nu_2 = coarse.theta{3}(12,2)\), and \(\nu_3 = coarse.theta{3}(12,3)\).
+
+
+
In addition, for the treatment of the boundary conditions, the boundary is divided into two sections. The first is \(\Gamma \backslash \Gamma_{root}\), where Robin boundary conditions are applied; the second is \(\Gamma_{root}\), where the incoming heat flux is applied. For each segment in these sections, the associated global nodes are provided:
For each of the sections \(i\), the index \(m\) varies in the range \(m \in \{1, ..., s_i\}\), where the \(s_i\) are the number of segments in section \(i\). As an example, to find the nodes \(\nu_1\) and \(\nu_2\) for segment 5 in the first section, we would use \(\nu_1 = coarse.theta{6}(5,1)\) and \(\nu_2 = coarse.theta{6}(5,2)\).
+
+
+
To plot the temperature distribution, plotsolution.m can be used. If \(z \equiv u_h\) is the vector with the computed temperature values for each of the nodes, then a contour plot of the temperature distribution can be obtained by:
+
+
+
+
plotsolution(FE_grid.coarse, z)
+
+
+
+
The first argument is the mesh used in the calculation of \(z\), and the second is the solution vector \(z\).
+
+
+
For the storage of the finite element matrices, use Python’s sparse matrix data structure. Also, for the solution of the resulting linear systems, use the default solution methods provided in Python, i.e.:
+
+
+
+
u=A\F
+
+
+
+
to solve for the FEM solution \(u\).
+
+
+
In python, we have
+
+
+
+
Environment
+
+
+
+
+
import numpy as np
+import scipy as sp
+from scipy.sparse import csc_matrix
+from scipy.sparse.linalg import spsolve
+from scipy.io import loadmat
+import matplotlib.tri as tri
+import matplotlib.pyplot as plt
+
+
+
+
+
Load and plot the triangulation
+
+
+
+
+
grids = loadmat('FE_grid.mat',simplify_cells=True)
+print(grids.keys())
+coarse_grid = grids['FE_grid']['coarse']
+# show the keys in the grid
+print(coarse_grid.keys())
+print("number of nodes:",coarse_grid['nodes'])
+print("number of elements:",coarse_grid['elements'])
+x=coarse_grid['coor'][:,0]
+y=coarse_grid['coor'][:,1]
+z=np.sin(np.pi*x)*np.cos(np.pi*y)
+# be careful the indices must start at 0, in mat files they start at one, so substract 1
+triangles=np.concatenate(coarse_grid['theta'][0:5]-1)
+T=tri.Triangulation(x,y,triangles)
Setting up the development Environment
diff --git a/course-rom/quickstart.html b/course-rom/quickstart.html
index ed3404d..ccb910a 100644
--- a/course-rom/quickstart.html
+++ b/course-rom/quickstart.html
@@ -145,9 +145,10 @@
Setting up the development Environment
diff --git a/course-rom/tps/tp-assim-1.html b/course-rom/tps/tp-assim-1.html
index eb26e11..1eea2b0 100644
--- a/course-rom/tps/tp-assim-1.html
+++ b/course-rom/tps/tp-assim-1.html
@@ -145,9 +145,10 @@