pymor / pymor

pyMOR - Model Order Reduction with Python
https://pymor.org
Other
306 stars 106 forks source link

BTReductor.reduce with `projection='bfsr'` may return inconsistent rom #1971

Closed drittelhacker closed 1 year ago

drittelhacker commented 1 year ago

Describe the bug When run for an LTImodel for the CDPlayer benchmark ( E is the identity ) a rom with a non-identity E is genrated, but E is reported as E=NumpyMatrixOperator(<109x109 dense>, name='IdentityOperator') using projection='sr' or projection='biorth` avoids the problem

To Reproduce load the CDPlayer benchmark, generate LTIModel from matrices, generate the BTreductor and call reduce()

Expected behavior E=NumpyMatrixOperator(<109x109 dense>) is set and reported rather than indicating an identity that is not there.

System information: pyMOR Version 2022.2.0

Python 3.10.6 on Linux-5.19.0-35-generic-x86_64-with-glibc2.35

External Packages

DEALII: missing DUNEGDT: 2022.5.3 FENICS: missing GL: 3.1.6 IPYTHON: 8.4.1 IPYWIDGETS: 7.7.3 MATPLOTLIB: 3.7.0 MESHIO: 5.3.4 MPI: 3.1.4 NGSOLVE: missing NUMPY: 1.24.2 PYMESS: present PYTEST: 7.1.2 PYTHREEJS: 2.4.2 QT: PyQt5 (Qt 5.15.2) QTOPENGL: present SCIKIT_FEM: missing SCIPY: 1.10.1 SCIPY_LSMR: present SLYCOT: 0.5.3 SPHINX: 5.1.1 TORCH: 1.13.1+cu117 TYPER: 0.7.0 VTKIO: present

Additional context Add any other context about the problem here.

pmli commented 1 year ago

I can reproduce this on a simple example:

import numpy as np
from pymor.models.iosys import LTIModel
from pymor.reductors.bt import BTReductor

A = np.diag([-1, -2])
B = np.ones((2, 1))
C = np.ones((1, 2))
fom = LTIModel.from_matrices(A, B, C)
bt = BTReductor(fom)
rom = bt.reduce(1)
print(repr(fom))
print(repr(rom))

which prints

00:01 LTIPGReductor: Operator projection ...
00:01 LTIPGReductor: Building ROM ...
LTIModel(
    NumpyMatrixOperator(<2x2 dense>, source_id='STATE', range_id='STATE'),
    NumpyMatrixOperator(<2x1 dense>, range_id='STATE'),
    NumpyMatrixOperator(<1x2 dense>, source_id='STATE'),
    D=ZeroOperator(NumpyVectorSpace(1), NumpyVectorSpace(1)),
    E=IdentityOperator(NumpyVectorSpace(2, id='STATE')),
    presets={})
LTIModel(
    NumpyMatrixOperator(<1x1 dense>),
    NumpyMatrixOperator(<1x1 dense>),
    NumpyMatrixOperator(<1x1 dense>),
    D=ZeroOperator(NumpyVectorSpace(1), NumpyVectorSpace(1)),
    E=NumpyMatrixOperator(<1x1 dense>, name='IdentityOperator'),
    presets={},
    name='LTIModel_reduced')
pmli commented 1 year ago

The issue seems to be in this line:

https://github.com/pymor/pymor/blob/0ea528e10b08e7fbd97645667b5213b6be9525bf/src/pymor/algorithms/projection.py#L147

Because op is the IdentityOperator, the name of the projected Operator becomes 'IdentityOperator'.

How about replacing name=op.name by name=op.name + '_projected'?

drittelhacker commented 1 year ago

Is this actually valuable information anywhere, or is there more potential harm than benefit?

pmli commented 1 year ago

@sdrave proposed removing name inheritance from the project method (adding _projected pollutes things). We might consider removing names at a later point.

I can work on a PR.