Importing Snapshots from OpenFOAM

This notebook uses pyforce to import some parametric data generated with OpenFOAM-6, exploiting the fluidfoam library.

The snapshots are related to the Buoyant Cavity problem in fluid dynamics, governed by the Navier-Stokes equations, including energy, under the Boussinesq approximation. In particular, the snapshots have been generated using the case reported in ROSE-ROM4FOAM tutorials.

import numpy as np
import pickle
import matplotlib.pyplot as plt

from mpi4py import MPI
from dolfinx.io import gmshio
import gmsh
from dolfinx.fem import Function, FunctionSpace
import ufl

from pyforce.tools.functions_list import FunctionsList
from pyforce.tools.backends import LoopProgress
from pyforce.tools.write_read import import_OF, StoreFunctionsList as store

path = './Snapshots/OpenFOAM/'

importing = ['TrainSet', 'TestSet']

var_names = ['p', 'T', 'U']
is_vector = [False, False, True]
# var_names = ['U']
# is_vector = [True]

The snapshots are dependent on two different parameters: the Reynolds and the Richardson number, split into train and test set \begin{equation} \begin{array}{cc} Re_{train} = [15:5:150] & Ri_{train} = [0.2:0.4:5] \\ Re_{test} = [17.5:10:147.5] & Ri_{test} = [0.4:0.8:44] \end{array} \end{equation}

dRe = 5.
dRi = 0.4

# Train/Test Parameters
Re_train_test = [np.arange(15,  150+dRe/2, dRe), np.arange(15+dRe/2,  150+dRe/2, dRe*2)]
Ri_train_test = [np.arange(0.2,   5+dRi/2, dRi), np.arange(0.2+dRi/2,   5+dRi/2, dRi*2)  ]

mu_train_test = [np.meshgrid(Re_train_test[kk], Ri_train_test[kk]) for kk in range(len(importing))]

fig = plt.figure(figsize=(6,6))

plt.plot(mu_train_test[0][0].flatten(), mu_train_test[0][1].flatten(), 'ro', label='Train')
plt.plot(mu_train_test[1][0].flatten(), mu_train_test[1][1].flatten(), 'gs', label='Test')
plt.xlabel('Reynolds Number $Re$', fontsize=20)
plt.ylabel('Richardson Number $Ri$', fontsize=20)
plt.grid()
plt.legend(framealpha=1, fontsize=15, loc='upper right')

<matplotlib.legend.Legend at 0x7fc379e6d7e0>

Let us generate the mesh for importing OpenFOAM dataset into dolfinx

mesh_comm = MPI.COMM_WORLD
model_rank = 0

# Initialize the gmsh module
gmsh.initialize()

# Load the .geo file
gmsh.merge('cavity.geo')
gmsh.model.geo.synchronize()

# Set algorithm (adaptive = 1, Frontal-Delaunay = 6)
gmsh.option.setNumber("Mesh.Algorithm", 6)
gdim = 2

# Linear Finite Element
gmsh.model.mesh.generate(gdim)
gmsh.model.mesh.optimize("Netgen")

# Import into dolfinx
model_rank = 0
domain, ct, ft = gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, model_rank, gdim = gdim )
gmsh.finalize()

########################################################################################################

tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)

Info    : Reading 'cavity.geo'...
Info    : Done reading 'cavity.geo'
Info    : Meshing 1D...
Info    : [  0%] Meshing curve 1 (Line)
Info    : [ 30%] Meshing curve 2 (Line)
Info    : [ 50%] Meshing curve 3 (Line)
Info    : [ 80%] Meshing curve 4 (Line)
Info    : Done meshing 1D (Wall 0.000607034s, CPU 0.000692s)
Info    : Meshing 2D...
Info    : Meshing surface 1 (Transfinite)
Info    : Done meshing 2D (Wall 0.00732976s, CPU 0.007374s)
Info    : 16384 nodes 32770 elements
Info    : Optimizing mesh (Netgen)...
Info    : Done optimizing mesh (Wall 3.26009e-06s, CPU 1e-05s)

Let us define the functional space onto which the OpenFOAM data are projected.

fun_spaces = [FunctionSpace(domain, ('Lagrange', 1)), FunctionSpace(domain, ('Lagrange', 1)), FunctionSpace(domain, ufl.VectorElement("CG", domain.ufl_cell(), 1))]
# fun_spaces = [FunctionSpace(domain, ufl.VectorElement("CG", domain.ufl_cell(), 1))]

Let us import the pressure \(p\), temperature \(T\) and velocity \(\mathbf{u}\) fields. The mapping between OpenFOAM and dolfinx is performed using \(N\)-dimensional interpolation implemented in scipy.

impo_i = 0

dolfinx_path = './Snapshots/'+importing[impo_i]+'_'
snaps = {var_names[field_i]: FunctionsList(fun_spaces[field_i]) for field_i in range(len(var_names))}

for field_i in range(len(var_names)):
    bar = LoopProgress('Import '+var_names[field_i]+' - '+importing[impo_i], final = len(Re_train_test[impo_i]) * len(Ri_train_test[impo_i]))

    caseI = 0
    for Re_i in range(len(Re_train_test[impo_i])):
        Re = Re_train_test[impo_i][Re_i]
        for Ri_i in range(len(Ri_train_test[impo_i])):
            Ri = Ri_train_test[impo_i][Ri_i]

            path_ = path+importing[impo_i]+'/'+f'Case_{caseI+0:03}_Re{Re:.2f}_Ri{Ri:.2f}'

            oF = import_OF(path=path_, extract_dofs=True)
            of_snaps = oF.import_field(var_names[field_i], vector=is_vector[field_i], verbose=False)[0]

            # Projection in dolfinx
            dolfinx_snap = oF.foam_to_dolfinx(fun_spaces[field_i], of_snaps, variables=['x', 'y'],  cut_value = oF.of_dofs[2,0])

            for uu in dolfinx_snap._list:
                snaps[var_names[field_i]].append(uu)

            bar.update(1)
            caseI += 1
    del bar

    store(domain, snaps[var_names[field_i]], var_names[field_i], dolfinx_path+var_names[field_i])

Import U - TrainSet: 364.000 / 364.00 - 1.324 s/it

Let us normalise the temperature field, using min-max to have it scaled between \((0,1)\).

from pyforce.tools.write_read import ImportH5

field_i = 1
field = var_names[field_i]

for imp in importing:
    # imp = importing[0]
    dolfinx_path = './Snapshots/'+imp+'_'

    T_snaps = ImportH5(fun_spaces[field_i], dolfinx_path+field, field)[0]

    if imp == 'TrainSet':
        _min = min([np.min(snap) for snap in T_snaps._list])
        _max = max([np.max(snap) for snap in T_snaps._list])

    T_norm_snaps = FunctionsList(T_snaps.fun_space)
    for snap in T_snaps._list:
        T_norm_snaps.append((snap - _min) / (_max - _min))

    store(domain, T_norm_snaps, 'norm_'+field, dolfinx_path+'norm_'+field)