I am implementing a linear interpolation between nodal values to compare FOM and ROM solutions.

Since I had already implemented a probe interpolation using the fenics __call__ method for CoefficientFunctions, I expected my interpolation to match these values.

To my surprise, they don't.
In fact, they do at the outflow (which is reassuring although trivial. Reassuring because the mass conservation calculation uses this value; trivial because it is a nodal value, so no interpolation takes place.), but they don't at the mid-point x=0.5.

I think there might be a problem with the moving mesh and the interpolation routine used inside fenics.
I am going to stick to my interpolation and use it to compare FOM and ROM solutions.

legend:

• dashed: expected value.
• cont.: interpolation.

• Review BDF-2 definition in the manuscript.
• Implement.
• Test:
  • Test for existing BDF-1 solution.
  • Test for BDF-2 solution, use conservation of mass as a validation.

See references:

• integration-seance3.pdf

The lifting operator is introduced as a linear combination of the boundary conditions.

Originally, this was the expression used:

f = fenics.Expression(
f"{right} * (x[0] / L) + {left} * (L - x[0]) / L",
degree=2,
L=L,
dLt_dt=dLt_dt,
t=t,
**mu,
)

This is a fragile definition, because it assumes the user will define right and left with parentheses.

The fixup is to include the parentheses by default.

f = fenics.Expression(
f"({right}) * (x[0] / L) + ({left}) * (L - x[0]) / L",
degree=2,
L=L,
dLt_dt=dLt_dt,
t=t,
**mu,
)

This showed up when I included a moving domain, because the RHS boundary condition is a sum of terms.
So far all the boundary conditions were made out of products (sighs...)

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

Reference: KikeM/msc-thesis#51

I closed KikeM/msc-thesis#19 in PR-#34 because I thought I had correctly implemented the MDEIM for the mass matrix too. There was a bug in the code, the ROM interpolation was being computed and then overwritten by the FOM assembly and projection.

def assemble_mass(self, mu, t):
        
       if self.mdeim_Mh:            
            MN = self.mdeim_Mh.interpolate(mu=mu, t=t, which=self.mdeim_Mh.ROM)        
        
        # Assemble FOM operator        
        Mh = self.fom.assemble_mass(mu, t)        
        MN = self.to_rom_bilinear(Mh)
        
        return MN 

Given an initial condition in the complete problem, the initial condition for the homogenous problem is obtained by subtracting the initial value of the Dirichlet lifting.

This should be done in fom/base.py, right before starting time integration.

For the piston problem, we will only use this for the nonlinear term.

I think there is a bug with the BDF-2 scheme.

I am not sure if I am setting correctly the coefficients multiplying the forcing terms that arise due to time discretization.

KikeM/msc-literature-review#41

• On which mesh is this computed? Fixed or moving?

Related to #46

Apparently, I should be able to change the mesh nodes by simply operating on the mesh coordinates.

If this is true, as long as I take into account the mesh orientation (right to left),
I can manually introduce the displacement field
(and hence not requiring the use of interpolation).

mesh = UnitSquareMesh(10, 10)
for x in mesh.coordinates():
    x[0] += 0.1*x[1]

From: https://fenicsproject.org/qa/11133/moving-mesh-to-specific-points/

More evidence in the same direction: https://fenicsproject.discourse.group/t/coordinate-update/799/13

I wonder what is going on there ...

It's true that there are three nested loops:

• Parameter space
  • Timestemps
    • Basis vectors

Perhaps I could combine the latter two and use multiprocessing.

That being said, it has more or less acceptable times (30 mins por 10 parameters for each ROM/DEIM/MDEIM) to get some results and pour them in a document.

For the fixed domain problem, the mass matrix is constant for all parameters and time.
Therefore, the interpolation should be exact, since there is only one basis vector: the matrix itself.
This means that the errors between the naive ROM with or without hyperreduction of the mass matrix should be the same.
However, when I fixed KikeM/msc-thesis#38 and launched again the tests, I realised the error was higher than what I expected.

The only difference was the mass matrix, whose interpolation is now correctly used.
But if the error was different, it meant that the mass matrix was different.

To have a better understanding, I parametrised the tests in tests/test_mdeim.py to assemble and interpolate the mass matrix too (those tests build the FOM operator reusing training parameters, and for the simple linear scenarios I am dealing with now it means that the FOM operator should be the same).

To my surprise, it is not being correctly built.

The reason: the Dirichlet entry is being selected as the interpolation node, but since the local assembly does not impose the boundary conditions, the FOM value corresponds to the integral of the boundary basis function.
The value of this integral is never seen in the FOM operator during the snapshot collection step, since it is overridden by the Dirichlet boundary conditions.
As a consequence, when the theta function is computed, an incorrect value is used in the RHS of the linear system.

This is useful for the sensitivity analysis.

Several bugs coupled in one:

• The rhs DEIM object was being added to the RomConstructor.FORCING block, instead of RomConstructor.RHS
• The DiscreteEmpiricalInterpolation.interpolate method did not make use of the which variable.
• The deim_rhs collateral basis was not projected.

Eventually I want to have tests that certify the code base and scripts that any user can pick up to use the ROM.

In a way, the tests are the ingredients for the scripts, but I don't to duplicate code.
Additionally, I don't want to check the output of the script, simply that it is running.
Therefore, I should implement something like this to have pytest check that the script is up and running.

I observed during the implementation of the DEIM algorithm (#22) that when I increase the number of nodes, the two boundary point where coming out as interpolation nodes (which is not necessarily a problem).

Yet, when I checked the basis vectors, I realised that the homogeneous boundary conditions are not being honoured by all basis vectors.

I went to check the boundary conditions of the vectors in the naive ROM for MFP-1 ... and not all of them satisfy the boundary values either!
I performed a visual check on the basis vectors, only to find out that most of them are smooth up to basis vector 11.

I think this might have to do with using the snapshots as they are, without any prior orthogonalisation, scaling with their norm, etc.

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.        ]])

Originally posted by @KikeM in https://github.com/KikeM/msc-thesis/issues/52#issuecomment-873440801

SVD honours constant values when they are zero, but not when they are 1.0.

Thus, the boundary conditions are not being honoured in the MDEIMs but they are in the DEIMs.

Originally posted by @KikeM in #14 (comment)

To compute the a posteriori error, we use another ROM with an additional basis elements.

The continuous L2 norm, which is an integral over the domain

translates into the discrete L2 norm divided by the number of elements in the vector (which is what dx represents),

The bug was that I was dividing by the number of elements in the vector, not the square root.