Projection-based ROM example

[dreamer2368]

Poisson Global ROM example

Demonstration result

Offline phase

python3 poisson_global_rom.py -offline -f 1.0 -id 0
python3 poisson_global_rom.py -offline -f 1.1 -id 1
python3 poisson_global_rom.py -offline -f 1.2 -id 2

Following is the output from one of the commands above. It takes longer than c++ libROM example, as we're using serial MFEM (details explained later)

Elapsed time for assembling FOM: 1.085698e+00 second

Elapsed time for solving FOM: 1.092814e+00 second

Merge phase

python3 poisson_global_rom.py -merge -ns 3

Following is the output from the merge command:

Opening file: basis0_snapshot.000000
Opening file: basis1_snapshot.000000
Opening file: basis2_snapshot.000000
Creating file: basis.000000
Elapsed time for merging and building ROM basis: 8.647323e-02 second

FOM phase

python3 poisson_global_rom.py -fom -f 1.15

The resulting output is similar to offline phase.

Online phase

python3 poisson_global_rom.py -online -f 1.15

Following is the output of the online command above:

Relative error of ROM solution = 6.40114E-04
Unable to connect to GLVis server at localhost:19916
GLVis visualization disabled.
Elapsed time for assembling ROM: 1.127186e+00 second

Elapsed time for solving ROM: 2.360344e-05 second

Differences from c++ libROM

1. Due to the installation issue of PyMFEM, current example uses serial MFEM classes instead of parallel classes. Solving time for FOM increased significantly as a result of not using AMG preconditioner.
2. ComputeCtAB defined in python side
   While there is ComputeCtAB provided from c++ libROM, this example defines its own ComputeCtAB, which is essentially a python copy of the same function. There are mainly two reasons why it is done this way:
   • The original ComputeCtAB from c++ libROM only supports parallel MFEM, which is currently not available due to the installation issue.
   • Even with a serial version implemented from c++ libROM, it is currently not possible to use this function, due to the incompatibility between SWIG and pybind11. While SWIG Object of type 'double *' has a workaround since double * is a fundamental data type, other objects such as mfem::HypreParMatrix or mfem::Operator do not have a way to convert from SWIG to pybind11.

Installation process updates

MFEM dependencies removed

• Originally, class/functions in libROM/lib/mfem require explicit access to MFEM class/functions. To bind them in pybind11, pylibROM now has explicit include/link to the mfem dependencies.
• However, due to the SWIG/pybind11 compatibility issue described above, even when these functions are properly binded, they cannot even be used in python, since we don't have access to the pointer inside SWIG object.
• Currently, the best way is to re-write libROM/lib/mfem functions purely in python, with combination of PyMFEM objects and pylibROM. In this case, libROM needs not be compiled with mfem dependencies, which simplifies the installation a little bit.

Pure python functions added

As a result, we re-implement these mfem-related functions purely in python. Currently two python functions are added:

• pylibROM.mfem.ComputeCtAB: compute projected reduced matrix from mfem::HypreParMatrix or mfem::Operator.
• pylibROM.python_utils.swigdouble2numpyarray: return numpy array from SWIG Object of type 'double *'

Previous __init__.py files were not written properly, and in fact not doing any job. Now they are properly re-written.

NOTE: usual cmake installation only includes c++ bindings from c++ libROM library, and thus these python modules are not available for cmake installation. They will be built into the module only with pip installation.

c++ bindings renamed as _pylibROM

In order to include both c++ bindings/python routines, the c++ bindings need to be renamed. Common practice is to prepand with an underscore: _pylibROM.

NOTE: usual cmake installation would have _pylibROM package, not pylibROM itself. If the package is compiled via cmake, then the users need to import _pylibROM instead of pylibROM.

Enforce memory contiguity of input numpy 1d array

Our current convention for double * is to handle them in the form of numpy.array. One issue with this is that numpy.array sometimes does not necessarily have a contiguous memory. libROM's implementation with double * is predicated upon the assumption that the array memory is always contiguous, thus requiring an enforcement at pybind interface.

For example, the numpy array col_vec in the following examples is contiguous in memory:
Valid case 1 Column vectors from fortran-contiguous 2d numpy array

import numpy as np
input_snapshot = np.array([[ 0.5377, -1.3077, -1.3499],
                           [ 1.8339, -0.4336,  3.0349],
                           [-2.2588,  0.3426,  0.7254],
                           [ 0.8622,  3.5784, -0.0631],
                           [ 0.3188,  2.7694,  0.7147]], order = 'F') # <- note the order = 'F'
idx = 1 # some random column index between 0 - 2
col_vec = input_snapshot[:, idx]

Valid case 2 The transpose of the matrix is constructed first then transposed back to the original shape. Its column vector is then contiguous in memory

input_snapshot = np.array([[0.5377, 1.8339, -2.2588, 0.8622, 0.3188],
                           [-1.3077, -0.4336, 0.3426, 3.5784, 2.7694],
                           [-1.3499, 3.0349, 0.7254, -0.0631, 0.7147]]) # <- note the order is default
input_snapshot = input_snapshot.T # <- note this transpose does not change the memory contiguity
idx = 1 # some random column index between 0 - 2
col_vec = input_snapshot[:, idx]

However, the following examples are not contiguous in memory:
Invalid case 1 Column vectors from the default 2d numpy array (row-major order)

input_snapshot = np.array([[ 0.5377, -1.3077, -1.3499],
                           [ 1.8339, -0.4336,  3.0349],
                           [-2.2588,  0.3426,  0.7254],
                           [ 0.8622,  3.5784, -0.0631],
                           [ 0.3188,  2.7694,  0.7147]]) # <- note the order is default
idx = 1 # some random column index between 0 - 2
col_vec = input_snapshot[:, idx]

Invalid case 2 Row vectors of the transposed matrix, which originally constructed in default row-major order

input_snapshot = np.array([[ 0.5377, -1.3077, -1.3499],
                           [ 1.8339, -0.4336,  3.0349],
                           [-2.2588,  0.3426,  0.7254],
                           [ 0.8622,  3.5784, -0.0631],
                           [ 0.3188,  2.7694,  0.7147]]) # <- note the order is default
input_snapshot = input_snapshot.T # <- note this transpose does not change the memory contiguity
idx = 1 # some random column index between 0 - 2
col_vec = input_snapshot[idx, :] # <- even though it's a row vector, it is still not contiguous in memory.

Now bindings/pylibROM/python_utils/cpp_utils.hpp provides auxiliary functions that ensure this contiguity. Mainly,

• double* getPointer(py::array_t<double> &u_in): returns double * from numpy array. If numpy array is 1D, then checks the contiguity.
• double* getVectorPointer(py::array_t<double> &u_in): returns double * from numpy array. Checks if the array is 1D and contiguous in memory.

All the libROM functions that take double * as an input 1d array are now simplified using the functions above. For example of DMD::takeSample,

.def("takeSample", [](DMD &self, py::array_t<double> &u_in, double t) {
    self.takeSample(getVectorPointer(u_in), t);
})

Minor fixes

• Previously, Database::formats was wrongly placed within BasisGenerator. This is now properly implemented with the pure class Database in utils.