diff --git a/README.md b/README.md index 61b4159..ffb06a5 100644 --- a/README.md +++ b/README.md @@ -4,7 +4,7 @@

-[![Reference Paper 1](https://img.shields.io/badge/Reference%20Paper%201-arXiv:%202401.07300-gray?labelColor=blue&style=flat&link=https://arxiv.org/abs/2401.07300)](https://arxiv.org/abs/2401.07300) [![Reference Paper 2](https://img.shields.io/badge/Reference%20Paper%202-10.1016/j.nucengdes.2024.113105-gray?labelColor=blue&style=flat&link=https://www.sciencedirect.com/science/article/pii/S002954932400205X)](https://www.sciencedirect.com/science/article/pii/S002954932400205X) [![Docs](https://img.shields.io/badge/Docs-green?style=flat&link=https://ermete-lab.github.io/ROSE-pyforce/intro.html)](https://ermete-lab.github.io/ROSE-pyforce/intro.html) [![Tutorials](https://img.shields.io/badge/Tutorials-red?style=flat&link=https://ermete-lab.github.io/ROSE-pyforce/tutorials.html)](https://ermete-lab.github.io/ROSE-pyforce/tutorials.html) +[![Reference Paper 1](https://img.shields.io/badge/Reference%20Paper%201-https://doi.org/10.1016/j.apm.2024.06.040-gray?labelColor=blue&style=flat&link=https://doi.org/10.1016/j.apm.2024.06.040)](https://doi.org/10.1016/j.apm.2024.06.040) [![Reference Paper 2](https://img.shields.io/badge/Reference%20Paper%202-10.1016/j.nucengdes.2024.113105-gray?labelColor=blue&style=flat&link=https://www.sciencedirect.com/science/article/pii/S002954932400205X)](https://www.sciencedirect.com/science/article/pii/S002954932400205X) [![Docs](https://img.shields.io/badge/Docs-green?style=flat&link=https://ermete-lab.github.io/ROSE-pyforce/intro.html)](https://ermete-lab.github.io/ROSE-pyforce/intro.html) [![Tutorials](https://img.shields.io/badge/Tutorials-red?style=flat&link=https://ermete-lab.github.io/ROSE-pyforce/tutorials.html)](https://ermete-lab.github.io/ROSE-pyforce/tutorials.html) [![Testing pyforce](https://github.com/Steriva/ROSE-pyforce/actions/workflows/testing.yml/badge.svg)](https://github.com/Steriva/ROSE-pyforce/actions/workflows/testing.yml) [![JOSS draft paper](https://github.com/ERMETE-Lab/ROSE-pyforce/actions/workflows/draft-pdf.yml/badge.svg)](https://github.com/ERMETE-Lab/ROSE-pyforce/actions/workflows/draft-pdf.yml) @@ -43,23 +43,27 @@ This package is aimed to be a valuable tool for other researchers, engineers, an If you are going to use *pyforce* in your research work, please cite the following articles. The authors would be pleased if you could cite the relevant papers: -1. Stefano Riva, Carolina Introini, and Antonio Cammi. Multi-physics model bias correction with data-driven reduced order modelling techniques: Application to nuclear case studies, 2024. [arXiv:2401.07300](http://arxiv.org/abs/2401.07300). +1. Stefano Riva, Carolina Introini, and Antonio Cammi, “Multi-physics model bias correction with data-driven reduced order techniques: Application to nuclear case studies,” Applied Mathematical Modelling, vol. 135, pp. 243–268, 2024. [https://doi.org/10.1016/j.apm.2024.06.040](https://doi.org/10.1016/j.apm.2024.06.040). 2. Antonio Cammi, Stefano Riva, Carolina Introini, Lorenzo Loi, and Enrico Padovani. Data-driven model order reduction for sensor positioning and indirect reconstruction with noisy data: Application to a circulating fuel reactor. Nuclear Engineering and Design, 421:113105, 2024. doi:[https://doi.org/10.1016/j.nucengdes.2024.113105](https://doi.org/10.1016/j.nucengdes.2024.113105). For LaTeX users: ```bibtex -@misc{RMP_2024, - title={Multi-Physics Model Bias Correction with Data-Driven Reduced Order Modelling Techniques: Application to Nuclear Case Studies}, - author={Stefano Riva and Carolina Introini and Antonio Cammi}, - year={2024}, - eprint={2401.07300}, - archivePrefix={arXiv}, - primaryClass={math.NA} +@article{RIVA2024_AMM, +title = {Multi-physics model bias correction with data-driven reduced order techniques: Application to nuclear case studies}, +journal = {Applied Mathematical Modelling}, +volume = {135}, +pages = {243-268}, +year = {2024}, +issn = {0307-904X}, +doi = {https://doi.org/10.1016/j.apm.2024.06.040}, +url = {https://www.sciencedirect.com/science/article/pii/S0307904X24003196}, +author = {Stefano Riva and Carolina Introini and Antonio Cammi}, +keywords = {Reduced order modelling, Data driven, Nuclear reactors, Multi-physics, Model correction}, } -@article{DDMOR_CFR, +@article{CAMMI2024_NED, title = {Data-driven model order reduction for sensor positioning and indirect reconstruction with noisy data: Application to a Circulating Fuel Reactor}, journal = {Nuclear Engineering and Design}, volume = {421}, diff --git a/docs/intro.md b/docs/intro.md index d06359e..88d915e 100644 --- a/docs/intro.md +++ b/docs/intro.md @@ -4,7 +4,7 @@ ![pyforce](images/immy_pyforce2.png) -[![Reference Paper 1](https://img.shields.io/badge/Reference%20Paper%201-arXiv:%202401.07300-gray?labelColor=blue&style=flat&link=https://arxiv.org/abs/2401.07300)](https://arxiv.org/abs/2401.07300) [![Reference Paper 2](https://img.shields.io/badge/Reference%20Paper%202-10.1016/j.nucengdes.2024.113105-gray?labelColor=blue&style=flat&link=https://www.sciencedirect.com/science/article/pii/S002954932400205X)](https://www.sciencedirect.com/science/article/pii/S002954932400205X) +[![Reference Paper 1](https://img.shields.io/badge/Reference%20Paper%201-https://doi.org/10.1016/j.apm.2024.06.040-gray?labelColor=blue&style=flat&link=https://doi.org/10.1016/j.apm.2024.06.040)](https://doi.org/10.1016/j.apm.2024.06.040) [![Reference Paper 2](https://img.shields.io/badge/Reference%20Paper%202-10.1016/j.nucengdes.2024.113105-gray?labelColor=blue&style=flat&link=https://www.sciencedirect.com/science/article/pii/S002954932400205X)](https://www.sciencedirect.com/science/article/pii/S002954932400205X) ## Description @@ -26,26 +26,30 @@ This work has been carried out at the [Nuclear Reactors Group - ERMETE Lab](http ## How to cite pyforce -If you are going to use *pyforce* in your research work, please cite the following articles. +If you are going to use *pyforce* in your research work, please cite the following articles. The authors would be pleased if you could cite the relevant papers: -1. Stefano Riva, Carolina Introini, and Antonio Cammi. Multi-physics model bias correction with data-driven reduced order modelling techniques: Application to nuclear case studies, 2024. [arXiv:2401.07300](http://arxiv.org/abs/2401.07300). +1. Stefano Riva, Carolina Introini, and Antonio Cammi, “Multi-physics model bias correction with data-driven reduced order techniques: Application to nuclear case studies,” Applied Mathematical Modelling, vol. 135, pp. 243–268, 2024. [https://doi.org/10.1016/j.apm.2024.06.040](https://doi.org/10.1016/j.apm.2024.06.040). 2. Antonio Cammi, Stefano Riva, Carolina Introini, Lorenzo Loi, and Enrico Padovani. Data-driven model order reduction for sensor positioning and indirect reconstruction with noisy data: Application to a circulating fuel reactor. Nuclear Engineering and Design, 421:113105, 2024. doi:[https://doi.org/10.1016/j.nucengdes.2024.113105](https://doi.org/10.1016/j.nucengdes.2024.113105). For LaTeX users: ```bibtex -@misc{RMP_2024, - title={Multi-Physics Model Bias Correction with Data-Driven Reduced Order Modelling Techniques: Application to Nuclear Case Studies}, - author={Stefano Riva and Carolina Introini and Antonio Cammi}, - year={2024}, - eprint={2401.07300}, - archivePrefix={arXiv}, - primaryClass={math.NA} +@article{RIVA2024_AMM, +title = {Multi-physics model bias correction with data-driven reduced order techniques: Application to nuclear case studies}, +journal = {Applied Mathematical Modelling}, +volume = {135}, +pages = {243-268}, +year = {2024}, +issn = {0307-904X}, +doi = {https://doi.org/10.1016/j.apm.2024.06.040}, +url = {https://www.sciencedirect.com/science/article/pii/S0307904X24003196}, +author = {Stefano Riva and Carolina Introini and Antonio Cammi}, +keywords = {Reduced order modelling, Data driven, Nuclear reactors, Multi-physics, Model correction}, } -@article{DDMOR_CFR, +@article{CAMMI2024_NED, title = {Data-driven model order reduction for sensor positioning and indirect reconstruction with noisy data: Application to a Circulating Fuel Reactor}, journal = {Nuclear Engineering and Design}, volume = {421}, @@ -63,7 +67,7 @@ keywords = {Hybrid Data-Assimilation, Generalized Empirical Interpolation Method --- ## What is Reduced Order Modelling? -In scientific literature the expression Reduced Order Modelling (ROM) {cite:p}`Quarteroni2016, MadayChapter2020, Degen2020_conference` is related to a set of techniques devoted to the search for an optimal coordinate system onto which some parametric solutions of Partial Differential Equations (PDEs) -typically called High-Fidelity (HF) or Full Order (FOM) Model - can be represented. These methods are very useful in multi-query and real-time scenarios, when quick and efficient solutions of models are required, e.g. optimization, uncertainty quantification and inverse problems {cite:p}`Guo_Veroy2021, Degen2022`. Recently, with the developments in data-driven modelling a lot of interest in the combination of data and models have been raised. ROM offers new opportunities both to integrate the model with experimental data in real-time and to define methods of sensor positioning, by providing efficient tools to compress the prior knowledge about the system coming from the parametrized mathematical model into low-dimensional forms. +In scientific literature the expression Reduced Order Modelling (ROM) {cite:p}`Quarteroni2016, MadayChapter2020, Degen2020_conference` is related to a set of techniques devoted to the search for an optimal coordinate system onto which some parametric solutions of Partial Differential Equations (PDEs) -typically called High-Fidelity (HF) or Full Order (FOM) Model - can be represented. These methods are very useful in multi-query and real-time scenarios, when quick and efficient solutions of models are required, e.g. optimization, uncertainty quantification and inverse problems {cite:p}`Guo_Veroy2021, Degen2022`. Recently, with the developments in data-driven modelling a lot of interest in the combination of data and models have been raised. ROM offers new opportunities both to integrate the model with experimental data in real-time and to define methods of sensor positioning, by providing efficient tools to compress the prior knowledge about the system coming from the parametrized mathematical model into low-dimensional forms. ### Reduced Basis Methods Among all ROM methods, Reduced Basis (RB) {cite:p}`Quarteroni2014, Hesthaven2016, Degen2020_certifiedRB` are a well-established and widely used class of ROM techniques, which are based on an offline-online paradigm. In the offline stage, a set of RB functions {math}`\phi_n(\mathbf{x})` are derived from an ensemble of high-fidelity solutions, called *snapshots*, yielding a low dimensional space that retains the main features of the full-order model. Different approaches can be used to construct the reduced basis, such as the **greedy** algorithms {cite:p}`Maday2006` and the **POD** {cite:p}`Berkooz1993`. Regardless of the construction strategy, an approximation of the high-fidelity solution is sought during the online stage as a linear combination of the RB functions, i.e. @@ -88,4 +92,4 @@ The general structure of the DDROM methods is reported in the figure below, the ![pyforce](images/DA_ROM.png) The Full Order (FOM) solution is encoded into a reduced coordinate system {cite:p}`NENE2023_RMP, RMP_2024`, ensuring a dimensionality reduction by searching for the dominant physics into the system. This process of encoding can be done through POD, GEIM or other more advanced algorithms and in general represents the heaviest part from the computational point of view. Then, the reduced coordinate system is used to select the most important locations in the domain to place experimental sensors. -Once this is the Data Assimilation process can begin, by combining real measures and the background knowledge in the form of a reduced coordinate system so that an improved system of coordinate is found enabling us to go back to the FOM representation with a decoding process. \ No newline at end of file +Once this is the Data Assimilation process can begin, by combining real measures and the background knowledge in the form of a reduced coordinate system so that an improved system of coordinate is found enabling us to go back to the FOM representation with a decoding process. diff --git a/pyforce/pyforce/offline/geim.py b/pyforce/pyforce/offline/geim.py index f592627..d24351e 100644 --- a/pyforce/pyforce/offline/geim.py +++ b/pyforce/pyforce/offline/geim.py @@ -53,7 +53,9 @@ def offline(self, train_snap: FunctionsList, Mmax: int, User-defined available positions for the sensors, if `None` the positions are taken from the mesh elements. sampleEvery : int, optional (default = 5) If `xm` is not `None`, sampling rate for the selection of points from the mesh - + verbose : boolean, optional (Default = False) + If `True`, print of the progress is enabled. + Returns ---------- maxAbsErr : np.ndarray diff --git a/pyforce/pyforce/offline/pod.py b/pyforce/pyforce/offline/pod.py index 0ab8ba4..b9fd8fb 100644 --- a/pyforce/pyforce/offline/pod.py +++ b/pyforce/pyforce/offline/pod.py @@ -1,6 +1,6 @@ # Offline Phase: Proper Orthogonal Decomposition # Author: Stefano Riva, PhD Student, NRG, Politecnico di Milano -# Latest Code Update: 24 May 2024 +# Latest Code Update: 12 July 2024 # Latest Doc Update: 24 May 2024 import numpy as np @@ -374,10 +374,13 @@ def projection(self, snap: np.ndarray, N: int = None): if N is None: N = self.Nmax - # Vh_star = scipy.linalg.solve(np.diag(self.sing_vals), np.dot(self.modes.return_matrix().T, snap), - # assume_a='sym') - Vh_star = np.dot(np.linalg.inv(np.diag(self.sing_vals)), - np.dot(self.modes.return_matrix().T, snap)) + if isinstance(snap, FunctionsList): + _snap = snap.return_matrix() + else: + _snap = snap + + Vh_star = np.dot(np.linalg.inv(np.diag(self.sing_vals[:N])), + np.dot(self.modes.return_matrix().T[:N], _snap)) return Vh_star diff --git a/pyforce/pyforce/online/geim_synthetic.py b/pyforce/pyforce/online/geim.py similarity index 81% rename from pyforce/pyforce/online/geim_synthetic.py rename to pyforce/pyforce/online/geim.py index 558667d..3fd6e74 100644 --- a/pyforce/pyforce/online/geim_synthetic.py +++ b/pyforce/pyforce/online/geim.py @@ -1,4 +1,4 @@ -# Synthetic Online Phase: Generalised Empirical Interpolation Method +# Online Phase: Generalised Empirical Interpolation Method # Author: Stefano Riva, PhD Student, NRG, Politecnico di Milano # Latest Code Update: 05 March 2024 # Latest Doc Update: 05 March 2024 @@ -11,15 +11,15 @@ from pyforce.tools.functions_list import FunctionsList from pyforce.tools.timer import Timer -# GEIM online (synthetic measures) +# GEIM online (synthetic and real measures) class GEIM(): r""" - This class can be used to perform the online phase of the GEIM algorihtm for synthetic measures :math:`\mathbf{y}\in\mathbb{R}^M` obtained as evaluations of the magic sensors on the snapshot :math:`u(\mathbf{x};\,\boldsymbol{\mu})` as + This class can be used to perform the online phase of the GEIM algorihtm for synthetic and real measures :math:`\mathbf{y}\in\mathbb{R}^M` either obtained as evaluations of the magic sensors on the snapshot :math:`u(\mathbf{x};\,\boldsymbol{\mu})` as .. math:: y_m = v_m(u)+\varepsilon_m \qquad \qquad m = 1, \dots, M - given :math:`\varepsilon_m` random noise (either present or not). + given :math:`\varepsilon_m` random noise (either present or not), or by real experimental data on the physical system. Parameters ---------- @@ -269,4 +269,59 @@ def reconstruct(self, snap: np.ndarray, M, noise_value = None): resid = np.abs(snap - interp) - return interp, resid, computational_time \ No newline at end of file + return interp, resid, computational_time + + def real_reconstruct(self, measure: np.ndarray): + r""" + The interpolant given the `measure` vector :math:`\mathbf{y}` input is computed, by solving the GEIM linear system + + .. math:: + \mathbb{B}\boldsymbol{\beta} = \mathbf{y} + + then the interpolant is computed and returned + + .. math:: + \mathcal{I}_M(\mathbf{x}) = \sum_{m=1}^M \beta_m[u] \cdot q_m(\mathbf{x}) + + Parameters + ---------- + measure : np.ndarray + Measurement vector, shaped as :math:`M \times N_s`, given :math:`M` the number of sensors used and :math:`N_s` the number of parametric realisation. + + Returns + ---------- + interp : np.ndarray + Interpolant Field :math:`\mathcal{I}_M` of GEIM + computational_time : dict + Dictionary with the CPU time of the most relevant operations during the online phase. + + """ + + M, Ns = measure.shape + + if M > self.Mmax: + print('The maximum number of measures must not be higher than '+str(self.Mmax)+' --> set equal to '+str(self.Mmax)) + M = self.Mmax + + # Variables to store the computational times + computational_time = dict() + timing = Timer() + + interps = FunctionsList(self.V) + + for mu in range(Ns): + + y = measure[:, mu] + + # Solving the linear system + timing.start() + coeff = la.solve(self.B[:M, :M], y, lower = True) + computational_time['LinearSystem'] = timing.stop() + + # Compute the interpolant and residual + timing.start() + interps.append(self.mf.lin_combine(coeff)) + + computational_time['Reconstruction'] = timing.stop() + + return interps, computational_time \ No newline at end of file diff --git a/pyforce/pyforce/online/tr_geim_synthetic.py b/pyforce/pyforce/online/tr_geim.py similarity index 85% rename from pyforce/pyforce/online/tr_geim_synthetic.py rename to pyforce/pyforce/online/tr_geim.py index 8c555e0..421930b 100644 --- a/pyforce/pyforce/online/tr_geim_synthetic.py +++ b/pyforce/pyforce/online/tr_geim.py @@ -1,4 +1,4 @@ -# Synthetic Online Phase: Tikhonov-Regularisation Generalised Empirical Interpolation Method +# Online Phase: Tikhonov-Regularisation Generalised Empirical Interpolation Method # Author: Stefano Riva, PhD Student, NRG, Politecnico di Milano # Latest Code Update: 06 March 2024 # Latest Doc Update: 06 March 2024 @@ -13,15 +13,15 @@ from typing import Tuple -# TR-GEIM online (synthetic measures) +# TR-GEIM online (synthetic and real measures) class TRGEIM(): r""" - This class can be used to perform the online phase of the TR-GEIM algorihtm for synthetic measures :math:`\mathbf{y}\in\mathbb{R}^M` obtained as evaluations of the magic sensors on the snapshot :math:`u(\mathbf{x};\,\boldsymbol{\mu})` as + This class can be used to perform the online phase of the TR-GEIM algorihtm for synthetic and real measures :math:`\mathbf{y}\in\mathbb{R}^M` obtained as evaluations of the magic sensors on the snapshot :math:`u(\mathbf{x};\,\boldsymbol{\mu})` as .. math:: y_m = v_m(u)+\varepsilon_m \qquad \qquad m = 1, \dots, M - given :math:`\varepsilon_m` random noise. + given :math:`\varepsilon_m` random noise, or by real experimental data on the physical system. Parameters ---------- @@ -382,4 +382,64 @@ def hyperparameter_tuning(self, snaps: FunctionsList, noise_value: float, lambda_opt = lambda_star_samples[np.argmin(abs_err.mean(axis=0))] - return lambda_opt, lambda_star_samples, abs_err.mean(axis = 0) \ No newline at end of file + return lambda_opt, lambda_star_samples, abs_err.mean(axis = 0) + + def real_reconstruct(self, measure: np.ndarray, reg_param: float): + r""" + The interpolant given the `measure` vector :math:`\mathbf{y}` input is computed, by solving the GEIM linear system + + .. math:: + \left(\mathbb{B}^T\mathbb{B}+\lambda \mathbb{T}^T\mathbb{T}\right)\boldsymbol{\beta} = \mathbb{B}^T\mathbf{y}+\lambda \mathbb{T}^T\mathbb{T} \langle{\boldsymbol{\beta}}\rangle + + then the interpolant is computed and returned + + .. math:: + \mathcal{I}_M(\mathbf{x}) = \sum_{m=1}^M \beta_m[u] \cdot q_m(\mathbf{x}) + + Parameters + ---------- + measure : np.ndarray + Measurement vector, shaped as :math:`M \times N_s`, given :math:`M` the number of sensors used and :math:`N_s` the number of parametric realisation. + reg_param : float + Regularising parameter :math:`\lambda` + + Returns + ---------- + interp : np.ndarray + Interpolant Field :math:`\mathcal{I}_M` of GEIM + computational_time : dict + Dictionary with the CPU time of the most relevant operations during the online phase. + + """ + + M, Ns = measure.shape + + if M > self.Mmax: + print('The maximum number of measures must not be higher than '+str(self.Mmax)+' --> set equal to '+str(self.Mmax)) + M = self.Mmax + + # Variables to store the computational times + computational_time = dict() + timing = Timer() + + interps = FunctionsList(self.V) + + for mu in range(Ns): + + y = measure[:, mu] + + # Solving the linear system + timing.start() + sys_matrix = (self.B[:M, :M].T @ self.B[:M, :M]) + reg_param * (self.T[:M, :M].T @ self.T[:M, :M]) + rhs = np.dot(self.B[:M, :M].T, y) + reg_param * np.dot((self.T[:M, :M].T @ self.T[:M, :M]), self.mean_beta[:M]) + + coeff = la.solve(sys_matrix, rhs) + computational_time['LinearSystem'] = timing.stop() + + # Compute the interpolant and residual + timing.start() + interps.append(self.mf.lin_combine(coeff)) + + computational_time['Reconstruction'] = timing.stop() + + return interps, computational_time \ No newline at end of file diff --git a/pyforce/pyforce/tools/write_read.py b/pyforce/pyforce/tools/write_read.py index 4696768..f43b8ae 100644 --- a/pyforce/pyforce/tools/write_read.py +++ b/pyforce/pyforce/tools/write_read.py @@ -1,7 +1,7 @@ # I/O tools # Author: Stefano Riva, PhD Student, NRG, Politecnico di Milano -# Latest Code Update: 24 May 2024 -# Latest Doc Update: 25 June 2023 +# Latest Code Update: 05 July 2024 +# Latest Doc Update: 05 July 2024 # from pyforce.tools.functions_list import FunctionsList from .functions_list import FunctionsList @@ -13,6 +13,7 @@ import os import numpy as np import fluidfoam as of +import pandas as pd from scipy.interpolate import NearestNDInterpolator, LinearNDInterpolator def StoreFunctionsList(domain, snap: FunctionsList, var_name: str, filename: str, order = None): @@ -290,4 +291,4 @@ def foam_to_dolfinx(self, V: FunctionSpace, snaps: list, variables: list, if verbose is not None: progressBar.update(1, percentage=False) - return snap_dolfinx \ No newline at end of file + return snap_dolfinx diff --git a/pyforce/setup.py b/pyforce/setup.py index c0483bc..4158dad 100644 --- a/pyforce/setup.py +++ b/pyforce/setup.py @@ -7,7 +7,7 @@ setup( name='pyforce', - version='0.1.2', + version='0.1.3', description='Python Framework data-driven model Order Reduction for multi-physiCs problEms', long_description=README_MD, long_description_content_type="text/markdown",