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.
1.1. Part1 - Finite element approximation
We start by setting Feel++ environment. The results will be stored in the directory feelppdb
.
import feelpp
from course_rom import laplacian
from course_rom.forms import *
import numpy as np
import json
import os
d=os.getcwd()
print(f"directory={d}")
e=feelpp.Environment(['lap'],config=feelpp.localRepository("."))
Results
--------------------------------------------------------------------------- ModuleNotFoundError Traceback (most recent call last) File:2 1 import feelpp ----> 2 from course_rom import laplacian 3 from course_rom.forms import * 4 import numpy as np ModuleNotFoundError: No module named 'course_rom'
class Fin2D:
def __init__(self,level='coarse'):
self.level=level
feelpp.Environment.setConfigFile(f"{d}/src/cases/laplacian/fin/fin4/fin2d-{level}.cfg")
data = laplacian.loadSpecs(f"{d}/src/cases/laplacian/fin/fin2d-{level}.json")
self.lap=laplacian.get(dim=2,order=1)
self.lap.setSpecs(data)
self.lap.initialize()
print(f"loaded {level} mesh with {self.lap.mesh().numGlobalElements()} elements")
self.decompose()
self.a_thetas=[self.atheta0,self.atheta1,self.atheta2,self.atheta3,self.atheta4,self.atheta5]
self.f_thetas=[self.ftheta0]
def decompose(self):
self.Fq=[]
for i,mat in enumerate(['Gamma_root']):
self.Fq.append(self.lap.assembleFlux(markers=[mat],coeffs=1))
self.Aq=[]
for i,mat in enumerate(['Post', 'Fin_1', 'Fin_2', 'Fin_3', 'Fin_4']):
C=np.array([[1,0],[0,1]])
self.Aq.append(self.lap.assembleGradGrad(markers=[mat],coeffs=C))
for i,mat in enumerate(['Gamma_ext']):
self.Aq.append(self.lap.assembleMass(markers=[mat],coeffs=1))
def Qa(self):
return len(self.Aq)
def Qf(self):
return len(self.Fq)
def atheta0(self,mu): return 1
def atheta1(self,mu): return mu[0]
def atheta2(self,mu): return mu[1]
def atheta3(self,mu): return mu[2]
def atheta4(self,mu): return mu[3]
def atheta5(self,mu): return mu[4]
def ftheta0(self,mu): return 1
def a(self,mu):
return sum( (theta(mu) * aq for theta, aq in zip(self.a_thetas, self.Aq)),start=form2(test=self.lap.Xh(),trial=self.lap.Xh()))
def f(self,mu):
return sum( (theta(mu) * fq for theta, fq in zip(self.f_thetas, self.Fq)),start=form1(test=self.lap.Xh()))
def solve(self,mu):
u = self.lap.Xh().element()
fmu=self.f(mu)
self.a(mu).solve(solution=u,rhs=fmu,rebuild=True)
return u,fmu(u)
Results
We will need the following two sampling methods
def log_random_sample(bounds, num_samples=1):
# bounds is a list of tuples (a, b) for each dimension
log_bounds = [(np.log(a), np.log(b)) for a, b in bounds]
samples = [np.exp(np.random.uniform(log_a, log_b, num_samples)) for log_a, log_b in log_bounds]
return np.array(samples).T # Transpose to get samples in shape (num_samples, dimensions)
def log_equidistributed_sample(bounds, num_samples=1):
# bounds is a list of tuples (a, b) for each dimension
log_bounds = [(np.log(a), np.log(b)) for a, b in bounds]
samples = [np.exp(np.linspace(log_a,log_b,num_samples )) for log_a, log_b in log_bounds]
return np.array(samples).T # Transpose to get samples in shape (num_samples, dimensions)
Results
We will need the following grids
fin_coarse = Fin2D(level='coarse')
fin_medium = Fin2D(level='medium')
fin_fine = Fin2D(level='fine')
Results
--------------------------------------------------------------------------- AttributeError Traceback (most recent call last) File:1 ----> 1 fin_coarse = Fin2D(level='coarse') 2 fin_medium = Fin2D(level='medium') 3 fin_fine = Fin2D(level='fine') File :4, in Fin2D.__init__(self, level) 2 def __init__(self,level='coarse'): 3 self.level=level ----> 4 feelpp.Environment.setConfigFile(f"{d}/src/cases/laplacian/fin/fin4/fin2d-{level}.cfg") 5 data = laplacian.loadSpecs(f"{d}/src/cases/laplacian/fin/fin2d-{level}.json") 6 self.lap=laplacian.get(dim=2,order=1) AttributeError: module 'feelpp' has no attribute 'Environment'
We can print the size of the meshes using the following code.
for fin in [fin_coarse,fin_medium, fin_fine]:
print(f"[{fin.level}] nelts: {fin.lap.mesh().numGlobalElements()}")
sampling=log_equidistributed_sample([(1.5,1.5),(.01,1)], num_samples=20)
Troots=dict()
for fin in [fin_coarse,fin_medium, fin_fine]:
Troots[fin.level]=[]
for i,sample in enumerate(sampling):
mu=[sample[0],sample[0],sample[0],sample[0],sample[1]]
u,s=fin.solve(mu)
print(f"s_{fin.level}({mu})={s}")
Troots[fin.level].append(s)
print(f"sampling={sampling[:,1]}")
Results
--------------------------------------------------------------------------- NameError Traceback (most recent call last) File:1 ----> 1 for fin in [fin_coarse,fin_medium, fin_fine]: 2 print(f"[{fin.level}] nelts: {fin.lap.mesh().numGlobalElements()}") 4 sampling=log_equidistributed_sample([(1.5,1.5),(.01,1)], num_samples=20) NameError: name 'fin_coarse' is not defined
We can plot the results using plotly.
import plotly.graph_objects as go
import numpy as np
fig = go.Figure()
for fin in [fin_coarse,fin_medium, fin_fine]:
fig.add_trace(go.Scatter(
x=sampling[:,1], y=Troots[fin.level], mode='lines+markers', name=f'Mean_{fin.level}_Gamma_root'))
fig.update_layout(title='Mean Temperature', xaxis_title='mu', yaxis_title='Temperature')
fig.show()
Results
--------------------------------------------------------------------------- NameError Traceback (most recent call last) File:5 2 import numpy as np 4 fig = go.Figure() ----> 5 for fin in [fin_coarse,fin_medium, fin_fine]: 6 fig.add_trace(go.Scatter( 7 x=sampling[:,1], y=Troots[fin.level], mode='lines+markers', name=f'Mean_{fin.level}_Gamma_root')) 8 fig.update_layout(title='Mean Temperature', xaxis_title='mu', yaxis_title='Temperature') NameError: name 'fin_coarse' is not defined
1.2. Part 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,
and that the outputis given by
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
The goal is to implement an offline/ online version of the reduced – basis method following the computational decomposition indicated below.
-
Offline
-
Choose \(N\).
-
Choose the sample \(S_N\) .
-
Construct \(Z\).
-
Construct \(A^q_N, q = 1,\ldots,Q; F_N; \text{ and } L_N.\)
-
-
Online
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 commandsemilogy
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 commandsemilogy
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\).