Posts

THE POD-DL-ROM Script by Elise Grosjean (ENSTA-Paris, elise.grosjean@ensta-paris.fr) import sys #!{sys.executable} -m pip install numpy #!{sys.executable} -m pip install matplotlib #!{sys.executable} -m pip install scikit-fem #!{sys.executable} -m pip install tensorflow Let us consider a parameterized problem. Having (previously computed) solutions for different parameter values, we aim at approaching much faster “online” the solution associated to a new parameter value. To do so, we will employ a POD-DL-ROM strategy as in the paper “POD-DL-ROM: Enhancing deep learning-based reduced order models for nonlinear parametrized PDEs by proper orthogonal decomposition”. The idea is to first compress the training solutions with a POD before relying on a DL-ROM approach (see also “A comprehensive deep learning-based approach to reduced order modeling of nonlinear time-dependent parametrized PDEs”).

Dec 13, 2024

Let us consider a parameterized problem $\mathcal{P}$ we want to solve for different parameter values in a set $\mathcal{G}$. In general, Reduced Basis Methods (RBMs) are decomposed in two stages:

Dec 13, 2024

OFFLINE PART (SVD/POD/Greedy/rSVD…) Script by Elise Grosjean (elise.grosjean@ensta-paris.fr) import sys !{sys.executable} -m pip install numpy !{sys.executable} -m pip install matplotlib !{sys.executable} -m pip install scikit-fem [notice] A new release of pip is available: 24.1.1 -> 25.0.1 [notice] To update, run: python3.11 -m pip install --upgrade pip Let us consider a parameterized problem we want to solve for different parameter values in a set $\mathcal{G}$.

Dec 13, 2024

THE PBDW method import sys !{sys.executable} -m pip install numpy !{sys.executable} -m pip install matplotlib !{sys.executable} -m pip install scikit-fem Let us present one method that combines model order reduction and a data assimilation problem: the Parametrized-Background Data-Weak method (PBDW). $$ G^{bk,\mu}(u^{bk}(\mu)) = 0.$$$$ \mathcal{M}^{bk} = \{u^{bk}(\mu) : \mu \in \mathcal{P}^{bk}\}.$$The PBDW formulation integrates the parameterized mathematical model $G^{pb}$ and $M$ experimental observations associated with a parameter configuration $\mu^{true}$ to estimate the true field $u^{true}(\mu^{true})$ as well as any desired output $l_{out}(u^{true}(\mu^{true}\ )) \in \mathbb{R}$ for given output functional $l_{out}$. We intend that $\lVert u^{true}(\mu^{true})-u^{bk}(\mu^{true}) \rVert$ is small (i.e. that our model represents our data observations well).

Dec 13, 2024

import numpy as np import skfem # for Finite Element Method from skfem import Basis, ElementTriP0 from skfem import MeshTri import matplotlib.pyplot as plt import time # ============================================================ # Param # ============================================================ Nx = 50 Ny = 50 # Maillage triangulaire du carré (0,1)^2 m = MeshTri.init_tensor( np.linspace(0.0, 1.0, Nx + 1), np.linspace(0.0, 1.0, Ny + 1) ) dx=1./Nx dy=1./Ny cell_areas=0.5*dx*dy NumberOfModes = 15 NumberOfSnapshots = 20 MuMin = -2.0 MuMax = -1.0 rw = 0.15 # width of gaussian sensors basis = Basis(m, ElementTriP0()) # P0 : one value for each element dof_locs = basis.doflocs # barycenters of elements for P0 ndof = dof_locs.shape[1] # training set mu1_train = np.linspace(MuMin, MuMax, NumberOfSnapshots) mu2_train = np.linspace(MuMin, MuMax, NumberOfSnapshots) # ============================================================ # Function g(x; mu1, mu2) # ============================================================ def g_fun(x, mu1, mu2): return 1.0 / ((x[0] - mu1)**2 + (x[1] - mu2)) def project_g(mu1, mu2): return basis.project(lambda x: g_fun(x, mu1, mu2)) # ============================================================ # Gaussian kernel ell_m(x) # ============================================================ def gaussian_weight_noC(x, xm, rw): """ Gaussian centered at xm, without normalization constant. x: shape (2, npts) xm: shape (2,) """ return np.exp(-0.5 * np.sum((x - xm[:, None])**2, axis=0) / rw**2) # ============================================================ # Construction of continuous-inspired candidate sensors l_m # ============================================================ def sigma_m_from_center(xm, rw): """ Build the discrete representation of the continuous functional sigma_m(g) = \int_\Omega l_m(x) g dx for P0 functions, using barycenter quadrature: \int_K l_m(x) dx ~ |K| l_m(x_K)=w[K] Returns a weight vector w such that sigma_m(g_h) ~ sum_K w[K] * g_h[K] """ l_vals = gaussian_weight_noC(dof_locs, xm, rw) # l(x_K) w_K = cell_areas * l_vals # |K| * l(x_K) s = np.sum(w_K) if s > 1e-14: w_K = w_K / s # normalize so that l_m(1) = 1 approximately return w_K def apply_sigma_m(w, u): """ Apply the functional represented by weights w to the discrete field u. """ return np.dot(w, u) # candidate sensors: one candidate centered at each barycenter L_candidates = [] L_centers = [] for k in range(ndof): xm = dof_locs[:, k] wm = sigma_m_from_center(xm, rw) L_candidates.append(wm) L_centers.append((xm[0], xm[1])) # ============================================================ # OFFLINE : initialization # ============================================================ ns = NumberOfSnapshots ntrain = ns * ns GGrid = np.zeros((ndof, ntrain)) GinftyGrid = np.zeros(ntrain) mu_pairs = [] col = 0 for a in mu1_train: for b in mu2_train: u = project_g(a, b) GGrid[:, col] = u GinftyGrid[col] = np.linalg.norm(u, np.inf) mu_pairs.append((a, b)) col += 1 # first parameter muIndex = np.argmax(GinftyGrid) muIndexlist = [muIndex] mu1_star, mu2_star = mu_pairs[muIndex] g1 = GGrid[:, muIndex].copy() # first magic sensor l_1: maximize |l(g1)| scores = np.array([np.abs(apply_sigma_m(wm, g1)) for wm in L_candidates]) LIndex = np.argmax(scores) LIndexlist = [LIndex] Llist = [L_candidates[LIndex]] print("First magic sensor center:", L_centers[LIndex]) # first basis function q1 normalized by l_1(q1)=1 q1 = g1 / apply_sigma_m(Llist[0], g1) qlist = [q1] # ============================================================ # OFFLINE : GEIM iterations # ============================================================ for M in range(1, NumberOfModes): print(f"M = {M}") # Interpolation matrix B_M with entries l_i(q_j) Qmat = np.zeros((M, M)) for i in range(M): for j in range(M): Qmat[i, j] = apply_sigma_m(Llist[i], qlist[j]) ResidualGrid = np.zeros((ndof, ntrain)) ResidualInf = np.zeros(ntrain) # residual for each training parameter for col, (a, b) in enumerate(mu_pairs): u = GGrid[:, col].copy() rhs = np.array([apply_sigma_m(Llist[i], u) for i in range(M)]) alpha = np.linalg.solve(Qmat, rhs) residual = u.copy() for i in range(M): residual -= alpha[i] * qlist[i] ResidualGrid[:, col] = residual ResidualInf[col] = np.linalg.norm(residual, np.inf) # avoid already selected parameters for old_idx in muIndexlist: ResidualInf[old_idx] = -1e30 # new parameter muIndex = np.argmax(ResidualInf) muIndexlist.append(muIndex) mu1_star, mu2_star = mu_pairs[muIndex] rM = ResidualGrid[:, muIndex].copy() print("Selected parameter:", (mu1_star, mu2_star)) print("Residual max norm :", np.max(ResidualInf)) # choose new magic sensor maximizing |l(rM)| scores = np.array([abs(apply_sigma_m(w, rM)) for w in L_candidates]) # avoid already selected sensors for old_L in LIndexlist: scores[old_L] = -np.inf LIndex = np.argmax(scores) LIndexlist.append(LIndex) Llist.append(L_candidates[LIndex]) print("New magic sensor center:", L_centers[LIndex]) # normalize new basis function by l_M(q_M)=1 qM = rM / apply_sigma_m(Llist[-1], rM) qlist.append(qM) # ============================================================ # ONLINE # ============================================================ mu1_target = -0.5 mu2_target = -0.5 u_true = project_g(mu1_target, mu2_target) # complete interpolation matrix Qmat = np.zeros((NumberOfModes, NumberOfModes)) for i in range(NumberOfModes): for j in range(NumberOfModes): Qmat[i, j] = apply_sigma_m(Llist[i], qlist[j]) # rhs = [l_i(u_true)] rhs = np.array([apply_sigma_m(Llist[i], u_true) for i in range(NumberOfModes)]) # coefficients GEIM alpha = np.linalg.solve(Qmat, rhs) # reconstruction GEIM u_geim = np.zeros(ndof) for j in range(NumberOfModes): u_geim += alpha[j] * qlist[j] # ============================================================ # Test in one point # ============================================================ xseek, yseek = 0.5, 0.5 element_finder = m.element_finder() element_index = element_finder(np.array([xseek]), np.array([yseek]))[0] approx_value = u_geim[element_index] true_value_fem = u_true[element_index] true_value_exact = g_fun(np.array([[xseek], [yseek]]), mu1_target, mu2_target)[0] print() print(f"Tested point : ({xseek}, {yseek})") print(f"Approximation GEIM : {approx_value}") print(f"Projected FEM : {true_value_fem}") print(f"Exact value g(x,mu) : {true_value_exact}") print(f"Error |GEIM - FEM| : {abs(approx_value - true_value_fem)}") print(f"Error |GEIM - exact|: {abs(approx_value - true_value_exact)}") First magic sensor center: (0.01, 0.02) M = 1 Selected parameter: (-1.0, -2.0) Residual max norm : 0.02627306674745511 New magic sensor center: (1.0, 0.75) M = 2 Selected parameter: (-1.631578947368421, -1.0) Residual max norm : 0.015000046418788585 New magic sensor center: (0.01, 1.0) M = 3 Selected parameter: (-2.0, -2.0) Residual max norm : 0.0019474457678794344 New magic sensor center: (0.56, 0.01) M = 4 Selected parameter: (-1.2105263157894737, -1.0) Residual max norm : 0.00031219133964447207 New magic sensor center: (0.01, 0.46) M = 5 Selected parameter: (-1.0, -1.368421052631579) Residual max norm : 0.00026240204454349336 New magic sensor center: (1.0, 0.01) M = 6 Selected parameter: (-1.3157894736842106, -1.7894736842105263) Residual max norm : 5.765112454682722e-05 New magic sensor center: (0.43, 0.74) M = 7 Selected parameter: (-2.0, -1.105263157894737) Residual max norm : 1.5172047617061872e-05 New magic sensor center: (0.26, 0.01) M = 8 Selected parameter: (-1.0, -1.6842105263157894) Residual max norm : 6.4691190488781556e-06 New magic sensor center: (0.01, 0.7000000000000001) M = 9 Selected parameter: (-1.1578947368421053, -1.3157894736842106) Residual max norm : 1.9236475368019866e-06 New magic sensor center: (0.28, 0.43) M = 10 Selected parameter: (-1.1578947368421053, -2.0) Residual max norm : 1.5631815551238951e-06 New magic sensor center: (0.6900000000000001, 0.48) M = 11 Selected parameter: (-1.0, -1.105263157894737) Residual max norm : 3.9161614561813685e-07 New magic sensor center: (0.99, 1.0) M = 12 Selected parameter: (-1.105263157894737, -1.0) Residual max norm : 7.386577928874025e-07 New magic sensor center: (0.1, 0.21000000000000002) M = 13 Selected parameter: (-1.631578947368421, -2.0) Residual max norm : 1.7578529701250268e-07 New magic sensor center: (0.65, 1.0) M = 14 Selected parameter: (-2.0, -1.5789473684210527) Residual max norm : 1.6866644894449681e-07 New magic sensor center: (0.74, 0.01) Tested point : (0.5, 0.5) Approximation GEIM : 0.498339657743457 Projected FEM : 0.49834174836073053 Exact value g(x,mu) : 0.5 Error |GEIM - FEM| : 2.0906172735202233e-06 Error |GEIM - exact|: 0.0016603422565429904

Dec 13, 2024

THE POD-I METHOD Script by Elise Grosjean (ENSTA-Paris, elise.grosjean@ensta-paris.fr) import sys !{sys.executable} -m pip install numpy !{sys.executable} -m pip install matplotlib !{sys.executable} -m pip install scikit-fem !{sys.executable} -m pip install scikit-learn Let us consider a parameterized problem. Having (previously computed) solutions for different parameter values, we aim at approaching much faster the solution associated to a new parameter value. The Galerkin-Proper Orthogonal Decomposition (POD) algorithm is decomposed in two parts:

Dec 13, 2024

THE POD-GALERKIN METHOD Script by Elise Grosjean (ENSTA-Paris, elise.grosjean@ensta-paris.fr) import sys !{sys.executable} -m pip install numpy !{sys.executable} -m pip install matplotlib !{sys.executable} -m pip install scikit-fem Let us consider a parameterized problem. Having (previously computed) solutions for different parameter values, we aim at approaching much faster the solution associated to a new parameter value. The Galerkin-Proper Orthogonal Decomposition (POD) algorithm is decomposed in two parts:

Dec 13, 2024

THE PBDW method import sys !{sys.executable} -m pip install numpy !{sys.executable} -m pip install matplotlib !{sys.executable} -m pip install scikit-fem Let us present one method that combines model order reduction and a data assimilation problem: the Parametrized-Background Data-Weak method (PBDW). $$ G^{bk,\mu}(u^{bk}(\mu)) = 0.$$$$ \mathcal{M}^{bk} = \{u^{bk}(\mu) : \mu \in \mathcal{P}^{bk}\}.$$The PBDW formulation integrates the parameterized mathematical model $G^{pb}$ and $M$ experimental observations associated with a parameter configuration $\mu^{true}$ to estimate the true field $u^{true}(\mu^{true})$ as well as any desired output $l_{out}(u^{true}(\mu^{true}\ )) \in \mathbb{R}$ for given output functional $l_{out}$. We intend that $\lVert u^{true}(\mu^{true})-u^{bk}(\mu^{true}) \rVert$ is small (i.e. that our model represents our data observations well).

Dec 13, 2024

THE PBDW method import sys !{sys.executable} -m pip install numpy !{sys.executable} -m pip install matplotlib !{sys.executable} -m pip install scikit-fem Let us present one method that combines model order reduction and a data assimilation problem: the Parametrized-Background Data-Weak method (PBDW). $$ G^{bk,\mu}(u^{bk}(\mu)) = 0.$$$$ \mathcal{M}^{bk} = \{u^{bk}(\mu) : \mu \in \mathcal{P}^{bk}\}.$$The PBDW formulation integrates the parameterized mathematical model $G^{pb}$ and $M$ experimental observations associated with a parameter configuration $\mu^{true}$ to estimate the true field $u^{true}(\mu^{true})$ as well as any desired output $l_{out}(u^{true}(\mu^{true}\ )) \in \mathbb{R}$ for given output functional $l_{out}$. We intend that $\lVert u^{true}(\mu^{true})-u^{bk}(\mu^{true}) \rVert$ is small (i.e. that our model represents our data observations well).

Dec 13, 2024

THE NIRB TWO-GRID METHOD Script by Elise Grosjean (elise.grosjean@ensta-paris.fr) import sys !{sys.executable} -m pip install numpy !{sys.executable} -m pip install matplotlib !{sys.executable} -m pip install scikit-fem Let us consider a parameterized problem. The NIRB two-grid method follows the same ideas as in the Snapshot Proper Orthogonal Decomposition (POD) algorithm and is decomposed in two parts:

Dec 13, 2024