POD - compute correlation matrix

Stefano_Riva · October 10, 2022, 12:40pm

Hi everyone,
I’m trying to implement the Proper Orthogonal Decomposition in Fenicsx.
I have collected N_s snapshots u_i in a list. I need to compute the correlation matrix as

C_{ij} = \int_\Omega u_i\cdot u_j\,d\Omega

Then, the eigen-decomposition of this matrix must be performed C\vec{\eta}_n = \lambda_n\vec{\eta}_n to be used in the definition of the basis functions (POD modes) \phi_n

\phi_n(x) = \frac{1}{\sqrt{\lambda_n}}\sum_{k = 1}^{N_s}\eta_{n,k}\cdot u_n(x)

given \eta_{n,k} the k-th component of the n-th eigenvector.

In dolfinx, are there smart ways to do that? I had a look at the GitHub - RBniCS/RBniCSx at 92fc64f6a1493d2321c9e7f7dffe8ef7dadc0859, however I didn’t figure it out a simple way to compute the integral and to find the eigenvalues.

Any help would be helpful,
Thanks

francesco-ballarin · October 10, 2022, 3:58pm

In RBniCSx the relevant part related to the eigenvalue computation is at RBniCSx/proper_orthogonal_decomposition.py at main · RBniCS/RBniCSx · GitHub

Stefano_Riva · October 10, 2022, 6:59pm

Thank you very much for the quick reply.
The main issue is related to the definition of the correlation matrix. How is the compute_inner_product defined?

I had a look at the repository and to my understanding there are a lot of “dependencies”…

francesco-ballarin · October 10, 2022, 7:22pm

There are definitely dependencies, that’s just how a non trivial library works

You may want to look at the unit tests, which have a very basic and standalone example at

github.com

RBniCS/RBniCSx/blob/main/tests/unit/backends/test_backends_proper_orthogonal_decomposition.py#L84


      
              matrices = [dolfinx.fem.petsc.assemble_matrix(bilinear_form_cpp) for bilinear_form_cpp in bilinear_forms_cpp]
              for matrix in matrices:
                  matrix.assemble()
              tensors_list = rbnicsx.backends.TensorsList(bilinear_forms_cpp[0], mesh.comm)
              for matrix in matrices:
                  tensors_list.append(matrix)
              return tensors_list
          
          

          
@pytest.mark.parametrize("normalize", [True, False])
          def test_backends_proper_orthogonal_decomposition_functions(  # type: ignore[no-any-unimported]
              functions_list: rbnicsx.backends.FunctionsList, inner_product: ufl.Form, normalize: bool
          ) -> None:
              """Check rbnicsx.backends.proper_orthogonal_decomposition for the case of dolfinx.fem.Function snapshots."""
              compute_inner_product = rbnicsx.backends.bilinear_form_action(inner_product)
              eigenvalues, modes, eigenvectors = rbnicsx.backends.proper_orthogonal_decomposition(
                  functions_list[:2], compute_inner_product, N=2, tol=0.0, normalize=normalize)
              assert len(eigenvalues) == 2
              assert np.isclose(eigenvalues[0], 5)
              assert np.isclose(eigenvalues[1], 0)
              assert len(modes) == 2

From there, you will see that compute_inner_product is defined by passing a bilinear form to

github.com

RBniCS/RBniCSx/blob/main/rbnicsx/backends/projection.py#L493


      
                  If not provided, no postprocessing of the result will be carried out.
          
          
    Returns
              -------
              :
                  A list of callables that represents the action of L on a function.
              """
              return [linear_form_action(L_i) for L_i in L]
          
          

          
def bilinear_form_action(  # type: ignore[no-any-unimported]
              a: ufl.Form, part: typing.Optional[str] = None
          ) -> typing.Callable[[dolfinx.fem.Function], typing.Callable[[dolfinx.fem.Function], petsc4py.PETSc.ScalarType]]:
              """
              Return a callable that represents the action of a bilinear form on a pair of functions.
          
          
    Parameters
              ----------
              a
                  Bilinear form to be represented.
              part

and in particular

github.com

RBniCS/RBniCSx/blob/main/rbnicsx/backends/projection.py#L549


      
                          Function to be replaced to the test function.
          
          
            Returns
                      -------
                      :
                          Evaluation of the action of a on the provided pair of functions.
                      """
                      a_replacement_cpp.replace(test=fun_0)
                      return _extract_part(
                          a_replacement_cpp.comm.allreduce(
                              dolfinx.fem.assemble_scalar(a_replacement_cpp.form_cpp), op=mpi4py.MPI.SUM),
                          part)
          
          
        return _test_action
          
          
    return _trial_action
          
          

          
def block_diagonal_bilinear_form_action(  # type: ignore[no-any-unimported]
              a: typing.Sequence[ufl.Form], part: typing.Optional[str] = None
          ) -> typing.Sequence[

I would suggest that if you have further questions which are specific to RBniCSx rather than general FEniCSx, it would probably best to discuss them at

Topic		Replies	Views
Finding eigenvalues of a 3x3 tensor in dolfinx dolfinx	3	537	July 27, 2022
How to construct a sparse matrix used for calculating function values at multiple points dolfinx assembly , point-source	9	590	February 5, 2024
Copy of ufl.algebra.product	0	442	October 12, 2020
Fenicsx equivalent of compute_vertex_values() General	2	189	October 3, 2023
Corresponding functions in dolfinx	4	1789	April 17, 2024

POD - compute correlation matrix

Related Topics