ERMETE-Lab · Steriva · Jun 20, 2025 · Jun 16, 2025 · Jun 16, 2025 · Jun 20, 2025
diff --git a/.gitignore b/.gitignore
@@ -10,6 +10,10 @@ __pycache__/
 # C extensions
 *.so
 
+# Tutorial Data
+docs/Tutorials/04_FlowOverCylinder_OF2012/OF_case/
+docs/Tutorials/04_FlowOverCylinder_OF2012/Snapshots/
+
 # Distribution / packaging
 .Python
 build/

diff --git a/README.md b/README.md
@@ -15,7 +15,7 @@
 
 - [Description](#description)
 - [How to cite *pyforce*](#how-to-cite-pyforce)
-  - [Recent works with *pyforce*](#recent-works-with-pyforce)
+  - [Selected works with *pyforce*](#selected-works-with-pyforce)
 - [Installation](#installation)
 - [Tutorials](#tutorials)
 - [Authors and contributions](#authors-and-contributions)
@@ -119,7 +119,7 @@ The *pyforce* package is tested on some tutorials available in the [docs](https:
 
 *Coming Soon*: multiphysics (neutronics+thermal-hydraulics) with *dolfinx* and OpenFOAM.
 
-The snapshots can be either generated by the user or be downloaded at the following link [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.11483677.svg)](https://doi.org/10.5281/zenodo.11483677)
+The snapshots can be either generated by the user or be downloaded at the following link [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.13882851.svg)](https://zenodo.org/records/13882851)
 
 ## Authors and contributions
 

diff --git a/docs/Tutorials/02_MGDiffusion/03a_online_GEIM.ipynb b/docs/Tutorials/02_MGDiffusion/03a_online_GEIM.ipynb
@@ -226,12 +226,12 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 8,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
-    "from pyforce.online.geim_synthetic import GEIM\n",
-    "from pyforce.online.tr_geim_synthetic import TRGEIM\n",
+    "from pyforce.online.geim import GEIM\n",
+    "from pyforce.online.tr_geim import TRGEIM\n",
     "\n",
     "\n",
     "geim_online    = dict()\n",

diff --git a/docs/Tutorials/02_MGDiffusion/03b_online_PBDW.ipynb b/docs/Tutorials/02_MGDiffusion/03b_online_PBDW.ipynb
@@ -191,11 +191,11 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
-    "from pyforce.online.pbdw_synthetic import PBDW\n",
+    "from pyforce.online.pbdw import PBDW\n",
     "\n",
     "pbdw_online = dict()\n",
     "\n",

diff --git a/docs/Tutorials/03_BuoyantCavity_OF6/01_import_OFsnaps.ipynb b/docs/Tutorials/03_BuoyantCavity_OF6/01_import_OFsnaps.ipynb
diff --git a/docs/Tutorials/04_FlowOverCylinder_OF2012/01_import_snaps.ipynb b/docs/Tutorials/04_FlowOverCylinder_OF2012/01_import_snaps.ipynb
diff --git a/docs/Tutorials/04_FlowOverCylinder_OF2012/02_POD_GPR.ipynb b/docs/Tutorials/04_FlowOverCylinder_OF2012/02_POD_GPR.ipynb
diff --git a/docs/Tutorials/04_flow_cyl.rst b/docs/Tutorials/04_flow_cyl.rst
@@ -0,0 +1,26 @@
+Laminar Flow Past a Cylinder using OpenFOAM-v2012
+=============================================================
+
+**Aim of this tutorial:** learn how to import transient snapshots from OpenFOAM (version 2012 from com-version) using *pyvista* package, generate the POD modes from the snapshots for the Flow Over Cylinder problem using an accelerated version of the algorithm and create a reduced model using Gaussian Process Regression (GPR).
+
+The Flow Over Cylinder problem is a well-known benchmark in fluid dynamics, where a cylinder is placed in a flow field, leading to the formation of vortices and complex flow patterns. The problem is solved using OpenFOAM, a popular open-source CFD software (*pimpleFoam* solver in this case). The governing equations are the Navier Stokes equations:
+
+.. math::
+    \left\{
+        \begin{array}{ll}
+        \nabla\cdot \mathbf{u} = 0 & \mathbf{x}\in\Omega\\
+        \frac{\partial \mathbf{u}}{\partial t} + \left(\mathbf{u}\cdot \nabla\right)\mathbf{u} -\nu \Delta \mathbf{u} + \nabla p = 0 & \mathbf{x}\in\Omega\\
+        \mathbf{u}=\mathbf{u}_{in},\;\; \frac{\partial p}{\partial \mathbf{n}} = 0 & \mathbf{x}\in\Gamma_{in} \\
+        \mathbf{u}=\mathbf{0},\;\; \frac{\partial p}{\partial \mathbf{n}} = 0 & \mathbf{x}\in\Gamma_{w} \\
+        \frac{\partial \mathbf{u}}{\partial \mathbf{n}}=0,\;\; p = 0 & \mathbf{x}\in\Gamma_{out}
+        \end{array}
+    \right.
+
+given :math:`\Omega` as the domain and :math:`\partial\Omega` as its boundary, composed by :math:`\partial\Omega = \Gamma_{in}\cup\Gamma_{w}\cup\Gamma_{out}` where :math:`\Gamma_{in}` is the inlet boundary, :math:`\Gamma_{w}` is the wall boundary and :math:`\Gamma_{out}` is the outlet.
+
+.. toctree::
+    :maxdepth: 3
+    :caption: Steps:
+
+    Import Snapshots from OF <04_FlowOverCylinder_OF2012/01_import_snaps.ipynb>
+    Generation and testing of the surrogate mode: <04_FlowOverCylinder_OF2012/02_POD_GPR.ipynb>
diff --git a/docs/tutorials.rst b/docs/tutorials.rst
@@ -33,7 +33,11 @@ The following case studies are discussed:
     - *Offline*: Creation of the reduced space using the POD and GEIM.
     - *Online*:  Indirect Reconstruction (PE+POD-I), considering noisy data.
 
-*Coming Soon: multiphysics (neutronics+thermal-hydraulics) with dolfinx and OpenFOAM.*
+4. **Flow over Cylinder** (solved with OpenFOAM-v2012)
+
+    - *Offline*: Import from OpenFOAM and plotting using *pyvista*.
+    - *Offline*: Creation of the surrogate model for the dynamics of the reduced coefficients (POD) using Gaussian Process Regression.
+    - *Online*:  Indirect Reconstruction (PE+POD-I), considering noisy data.
 
 Each tutorial requires a dataset to be generated: the notebooks generating the dataset are available in the tutorials, alternatively the data can be downloaded from `Zenodo <https://doi.org/10.5281/zenodo.11483677>`_.
 
@@ -44,3 +48,4 @@ Each tutorial requires a dataset to be generated: the notebooks generating the d
     Unsteady Laminar Navier-Stokes - DFG2 benchmark <Tutorials/01_DFG2_benchmark.rst>
     MultiGroup Neutron Diffusion - ANL11-A2 benchmark <Tutorials/02_ANL11-A2_stationary.rst>
     Steady Buoyant Navier-Stokes - Differentially Heated Cavity <Tutorials/03_buoyant_cavity.rst>
+    Navier-Stokes (PimpleFoam) - Flow Over Cylinder <Tutorials/04_flow_cyl.rst>
diff --git a/pyforce/pyforce/offline/pod.py b/pyforce/pyforce/offline/pod.py
@@ -36,27 +36,59 @@ class POD():
       If `True`, print of the progress is enabled.
 
     """
-  def __init__(self, train_snap: FunctionsList, name: str, use_scipy=False, verbose = False) -> None:
+  def __init__(self, train_snap: FunctionsList, name: str, 
+               svd_acceleration_rank: int = None,
+               use_scipy=False, verbose = False) -> None:
 
     self.Ns   = len(train_snap)
     self.V = train_snap.fun_space
     self.norm = norms(self.V)
     self.name = name
 
     # Generate the correlation matrix using inner product in L^2
-    corrMatrix = np.zeros((self.Ns, self.Ns))
-    if verbose:
-        progressBar = LoopProgress(msg = "Computing " + self.name + ' correlation matrix', final = self.Ns)
+    if svd_acceleration_rank is not None:
+      _svd_u, _svd_s, _ = randomized_svd(train_snap.return_matrix(), n_components=svd_acceleration_rank, n_iter='auto')
 
-    for ii in range(self.Ns):
-        for jj in range(self.Ns):
-            if jj>=ii:
-                corrMatrix[ii,jj] = self.norm.L2innerProd(train_snap(ii), train_snap(jj))
-            else:
-                corrMatrix[ii,jj] = corrMatrix[jj,ii]
+      # Compute the residual energy and check if the rank is too low or too high
+      residual_energy = np.sum(_svd_s[:-1]**2) / np.sum(_svd_s**2)
+      if residual_energy <= 0.99:
+          print("Warning: The residual energy of the SVD is {} <= 0.9 for rank {}. This may indicate that the rank is too low.".format(residual_energy, svd_acceleration_rank))
+
+      _svd_v = _svd_u.T @ train_snap.return_matrix() # shape (svd_acceleration_rank, Ns)
 
-        if verbose:
-            progressBar.update(1, percentage = False)
+      # Compute the matrix for L2 inner product
+      _P_matrix = np.zeros((svd_acceleration_rank, svd_acceleration_rank))
+
+      if verbose:
+          progressBar = LoopProgress(msg = "Computing " + self.name + ' correlation matrix', final = svd_acceleration_rank)
+
+      for ii in range(svd_acceleration_rank):
+        for jj in range(svd_acceleration_rank):
+          if jj >= ii:
+            _P_matrix[ii, jj] = self.norm.L2innerProd(_svd_u[:, ii], _svd_u[:, jj])
+          else:
+            _P_matrix[ii, jj] = _P_matrix[jj, ii]
+        progressBar.update(1, percentage = False)
+
+      # Compute the correlation matrix using the SVD
+      corrMatrix = _svd_v.T @ _P_matrix @ _svd_v
+      assert corrMatrix.shape == (self.Ns, self.Ns), "The correlation matrix has the wrong shape: {}".format(corrMatrix.shape)
+
+    else:
+
+      if verbose:
+          progressBar = LoopProgress(msg = "Computing " + self.name + ' correlation matrix', final = self.Ns)
+
+      corrMatrix = np.zeros((self.Ns, self.Ns))
+      for ii in range(self.Ns):
+          for jj in range(self.Ns):
+              if jj>=ii:
+                  corrMatrix[ii,jj] = self.norm.L2innerProd(train_snap(ii), train_snap(jj))
+              else:
+                  corrMatrix[ii,jj] = corrMatrix[jj,ii]
+
+          if verbose:
+              progressBar.update(1, percentage = False)
 
     # Solving the eigenvalue problem and sorting the eigenvalue/eigenvector pairs
     if use_scipy:
@@ -67,9 +99,6 @@ def __init__(self, train_snap: FunctionsList, name: str, use_scipy=False, verbos
     eigenvalues = eigenvalues[sorted_indexes]
     eigenvectors = eigenvectors[:,sorted_indexes]
 
-    # if sum(eigenvalues < 0) > 0: 
-    #     warnings.warn("Check eigenvalues: some of them are negative! \n Try DiscretePOD with `svd` method")
-
     # Store the eigenvalue/eigenvector pairs
     self.eigenvalues  = eigenvalues
     self.eigenvectors = eigenvectors