BUG, MDEIM : Convection Operator

[KikeM]

What is going on with the HROM when the convection term is reduced too?

Reference: KikeM/msc-thesis#51

[KikeM]

I build the complete reduced operator by projection of the FOM operator and by addition of the approximation reduced operators.

To my surprise, they seem transposed!

(Pdb++) KN
array([[4.28326371, 0.01764899, 0.56785761],
       [0.05217947, 1.4404638 , 0.21454047],
       [0.52433937, 0.20503108, 7.47078799]])
(Pdb++) KN_m
array([[4.28326371, 0.05217947, 0.52433937],
       [0.01764899, 1.4404638 , 0.20503108],
       [0.56785761, 0.21454047, 7.47078799]])

[KikeM]

Ch = hrom.mdeim_convection.assemble(mu=mu_space[0], t=0.5)
Ch = bilinear_to_csr(Ch)
Ch = Ch.todense()

Ch_int = hrom.mdeim_convection.interpolate(
    mu=mu_space[0],
    t=0.5,
    which=OperatorType.FOM,
)
Ch_int = Ch_int.todense()

(Pdb++) Ch
matrix([[ 1.        ,  0.        ,  0.        ,  0.        ,  0.        ],
        [ 0.03140622, -0.00628124, -0.02512497,  0.        ,  0.        ],
        [ 0.        ,  0.02198435, -0.00628124, -0.01570311,  0.        ],
        [ 0.        ,  0.        ,  0.01256249, -0.00628124, -0.00628124],
        [ 0.        ,  0.        ,  0.        ,  0.        ,  1.        ]])
(Pdb++) Ch_int
matrix([[ 1.        ,  0.        ,  0.        ,  0.        ,  0.        ],
        [ 0.03140622, -0.00628124,  0.02198435,  0.        ,  0.        ],
        [ 0.        , -0.02512497, -0.00628124,  0.01256249,  0.        ],
        [ 0.        ,  0.        , -0.01570311, -0.00628124, -0.00628124],
        [ 0.        ,  0.        ,  0.        ,  0.        ,  1.        ]])

Aha!

The error can be traced back to the actual interpolation of the FOM.

What I do not understand then is how am I getting such a low online evaluation error ...

(Pdb++) hrom.mdeim_convection.errors_rom[0]
array([4.16546388e-17, 3.03180046e-17, 4.51616346e-17, 3.02895883e-17,
       6.38407379e-17, 4.03805237e-17, 5.04949207e-17, 4.03817555e-17,
       4.51427696e-17, 5.04834682e-17])

My guess is that probably in the vector form these two arrays are the same, but the bug is located when building back the matrix form. Let's have a look.

[KikeM]

for mu in tqdm(space, desc=msg_mu, leave=True):

            mu_idx, mu = self.add_mu(step=Stage.ONLINE, mu=mu)

            for t in tqdm(ts, desc=msg_time, leave=False):

                # Exact solution
                fh = self.assemble_snapshot(mu, t)

                #  Compute approximation
                fh_appr = self._interpolate(mu, t, which=self.FOM)

                breakpoint()

                error = self._compute_error(u=fh_appr, ue=fh)
                self.errors_rom[mu_idx].append(error)

            self.errors_rom[mu_idx] = np.array(self.errors_rom[mu_idx])

Indeed, in the vector form the two arrays are the same.

(Pdb++) fh
array([ 1.        ,  0.03164381, -0.00632876, -0.02531504,  0.02215066,
       -0.00632876, -0.0158219 ,  0.01265752, -0.00632876, -0.00632876,
        1.        ])
(Pdb++) fh_appr
array([ 1.        ,  0.03164381, -0.00632876, -0.02531504,  0.02215066,
       -0.00632876, -0.0158219 ,  0.01265752, -0.00632876, -0.00632876,
        1.        ])

[KikeM]

There is a bug in the determination of the rows and the columns.

def get_matrix_topology(self, mu, t):
    """Get mesh rows and columns arrangement.

    Parameters
    ----------
    mu : dict
    t : float

    Returns
    -------
    rows : np.array
    cols : np.array
    """

    Ah = self._assemble_matrix(mu, t)
    Ah = eliminate_zeros(Ah)
    Ah_dense = Ah.todense()
    rows, cols, _ = get_nonzero_entries(Ah)

    if self.name == OperatorType.CONVECTION:
        _list = list(zip(rows, cols, Ah.data))
        _list = sorted(_list, key=lambda x: x[0])

        _true_list = []
        for row in range(Ah_dense.shape[0]):
            for col in range(Ah_dense.shape[0]):
                value = Ah_dense[row, col]
                if np.isclose(value, 0.0):
                    continue
                else:
                    _true_list.append((row, col, value))

        breakpoint()

    return rows, cols

(Pdb++) pp _list
[(0, 0, 1.0),
 (1, 0, 0.03928432036567021),
 (1, 1, -0.011224091533048632),
 (1, 2, 0.02244818306609726), <---------------------------
 (2, 1, -0.028060228832621575),
 (2, 2, -0.011224091533048628),
 (2, 3, -0.011224091533048632),
 (3, 3, 1.0)]
(Pdb++) pp _true_list
[(0, 0, 1.0),
 (1, 0, 0.03928432036567021),
 (1, 1, -0.011224091533048632),
 (1, 2, -0.028060228832621575), <---------------------------
 (2, 1, 0.02244818306609726),
 (2, 2, -0.011224091533048628),
 (2, 3, -0.011224091533048632),
 (3, 3, 1.0)]
(Pdb++) Ah_dense
matrix([[ 1.        ,  0.        ,  0.        ,  0.        ],
        [ 0.03928432, -0.01122409, -0.02806023,  0.        ],
        [ 0.        ,  0.02244818, -0.01122409, -0.01122409],
        [ 0.        ,  0.        ,  0.        ,  1.        ]])