Instationary discretizer cannot handle non-homogeneous boundary and initial conditions properly

[SchleussJulia]

Describe the bug In case of "incompatible" initial and boundary values (see example below) the builtin discretizer/ the InstationaryModel fails to handle the boundary values correctly.

To Reproduce

import numpy as np
from pymor.basic import *

problem = InstationaryProblem(

    StationaryProblem(
        domain=RectDomain(top='dirichlet',bottom='dirichlet',left='dirichlet',right='dirichlet'),

        diffusion=ConstantFunction(1., dim_domain=2),

        rhs=ConstantFunction(0., dim_domain=2),

        dirichlet_data=ConstantFunction(0., dim_domain=2),

    ),

    T=0.1,

    initial_data=ConstantFunction(1., dim_domain=2),

)

n_h = 20
h = np.sqrt(2)/n_h
nt = 10

fom, data = discretize_instationary_cg(problem, diameter=h, nt=nt)

fom.visualize(fom.solve())

Expected behavior Zero boundary values would be expected. (In example, homogeneous initial conditions and dirichlet boundary values with constant value one do not work as well.)

Screenshots

System information: Tested on current pymor main.

[SchleussJulia]

Maybe the argument "dirichlet_clear_diag=True" in the assembly of the l2_0_product would do the job, see https://github.com/pymor/pymor/blob/aff70232b38996c94d009d6d9b6f6feef66de0dc/src/pymor/discretizers/builtin/cg.py#L1156

[sdrave]

That would make l2_0_product singular, which is bad when trying to compute Riesz representatives with it. We could do so for the mass operator. However, that would probably only work for implicit Euler. Another approach might be to always work in $H^1_0$. The caveat here would be that for time-dependent boundary data, one either would have to settle on the time-discretization already during spatial discretization in order to compute the correct RHS. Or one would have to also keep the full H^1-matrices (without boundary treatment) in memory to compute the right RHS during time-stepping. Any opinions @HenKlei, @ftschindler?