Reduced Jacobian of Linear Operator as (Fenics) Integral

[Frostyant]

I think I made a basic error here, bu if I have a reduced basis composed of exactly 1 vector and a linear problem, such as

F = (
            -df.inner(p,df.div(v))*df.dx +
            df.inner(df.grad(u),df.grad(v))*df.dx + 
            kappa*df.inner(u,v)*df.dx + 
            df.inner(q,df.div(u))*df.dx 
        )

Then the 1x1 reduced jacobian should just be the above integral but with u,p and v,q replaced by the reduced basis function.
Yet in the below code,

#!/usr/bin/env python
# This file is part of the pyMOR project (https://www.pymor.org).
# Copyright pyMOR developers and contributors. All rights reserved.
# License: BSD 2-Clause License (https://opensource.org/licenses/BSD-2-Clause)

from typer import Argument
from fenics import *
from sources.util.conversion import pymor_fenics

def alpha_mu(model,alpha):
    """
    from an alpha fenics function get the corresponding mu object and input it into fom
    """
    return  model.parameters.parse(list(alpha.vector()[:]))#translate mu as pymor parameter

def main(
    dim: int = Argument(..., help='Spatial dimension of the problem.'),
    n: int = Argument(..., help='Number of mesh intervals per spatial dimension.'),
    test: int = Argument(..., help='If True use TopologyFenicsOperator, else use FenicsOperator'),
):
    """Reduces a FEniCS-based nonlinear diffusion problem using POD/DEIM."""
    from pymor.tools import mpi

    
    fom,K,W,bcs,bc0s = discretize(dim, n, test)
    

    # ### ROM generation (POD/DEIM)
    from pymor.algorithms.newton import newton
    from pymor.algorithms.pod import pod
    from pymor.reductors.basic import StationaryRBReductor

    alphas = [1]
    U = fom.solution_space.empty()
    residuals = fom.solution_space.empty()
    for alpha_ in alphas:
        alpha = interpolate(Constant(alpha_),K)
        UU, data = newton(fom.operator, fom.rhs.as_vector(), mu=alpha_mu(fom,alpha), rtol=1e-6, return_residuals=True)
        U.append(UU)
        residuals.append(data['residuals'])

    alpha = interpolate(Constant(1.0),K)
    mu = alpha_mu(fom,alpha)

    rb, svals = pod(U, rtol=1e-7)
    reductor = StationaryRBReductor(fom, rb)
    rom = reductor.reduce()

    # ### ROM validation
    import time
    import numpy as np
    U = rom.operator.source.zeros()
    matrix = rom.operator.jacobian(U,mu=mu).matrix
    
    kappamin = Constant(1)
    kappamax = Constant(10)
    penalty = Constant(5)
    #Approximate alpha

    #Define kappa
    kappa = kappamax + (kappamin-kappamax)*alpha*(1+penalty)/(alpha+penalty)

    state_bases = rb
    matrix_2 = np.zeros([len(state_bases),len(state_bases)])
    for i, state_basis_1 in enumerate(state_bases):
        for j, state_basis_2 in enumerate(state_bases):
            vq = pymor_fenics(state_basis_2, W)
            vq_vector = vq.vector()
            up = pymor_fenics(state_basis_1, W)
            up_vector = up.vector()
            # Apply boundary conditions for state_basis_1
            for bc in bcs:
                bc.apply(up.vector())
            # Apply boundary conditions for state_basis_2
            for bc in bc0s:
                bc.apply(vq.vector())
            up = Function(W, up_vector)
            u, p = split(up)
            vq = Function(W, vq_vector)
            v, q = split(vq)
            linear_F = (
                    -inner(p, div(v)) * dx +  # Pressure Gradient term
                    inner(grad(u), grad(v)) * dx +  # Diffusion Term
                    kappa * inner(u, v) * dx +  # Topology term
                    inner(q, div(u)) * dx  # Divergence-Free Condition Term
            )

            matrix_2[i, j] = assemble(linear_F)
    print(matrix)
    print(matrix_2)


def discretize(dim, n, test=True):
    # ### problem definition
    import dolfin as df

    if dim == 2:
        mesh = df.UnitSquareMesh(n, n)
    else:
        raise NotImplementedError

    # FEM Spaces #
    FEMSpaces=[["Lagrange",2],["Lagrange",1],["Lagrange",1]]
    V = df.VectorElement(FEMSpaces[0][0],mesh.ufl_cell(),FEMSpaces[0][1])#velocity space
    Q = df.FiniteElement(FEMSpaces[1][0],mesh.ufl_cell(),FEMSpaces[1][1])#pressure space
    K = df.FunctionSpace(mesh,FEMSpaces[2][0],FEMSpaces[2][1])#function space of piecewise constants
    Welem = df.MixedElement([V,Q])
    W = df.FunctionSpace(mesh,Welem)
    Q = W.sub(1).collapse()
    # #

    def inflow_boundary(x, on_boundary):
        return near(x[0], 0) and (x[1] > 3.5/5) and (x[1] < 4.5/5)

    def outflow_boundary(x, on_boundary):
        return near(x[1], 0) and (x[0] > 3.5/5) and (x[0] < 4.5/5)

    def walls_boundary(x, on_boundary):
        return (near(x[0], 1) or near(x[1], 1) or
                (near(x[1], 0) and (x[0] <= 3.5/5 or x[0] >= 4.5/5)) or
                (near(x[0], 0) and (x[1] <= 3.5/5 or x[1] >= 4.5/5)))

    def DirichletBoundaries():
        """
        Defines boundaries
        """
        # Define boundary condition
        bcu_noslip = df.DirichletBC(W.sub(0), df.Constant((0, 0)), walls_boundary)
        bcu_inflow = df.DirichletBC(W.sub(0), df.Expression(("-(x[1]-3.5/5)*(x[1]-4.5/5)*(100)","0"), degree=2), inflow_boundary)  # parabolic input
        bcu_outflow = df.DirichletBC(W.sub(0), df.Expression(("0","-(x[0]-3.5/5)*(x[0]-4.5/5)*(-100)"), degree=2), outflow_boundary)  # parabolic input

        bdc = [bcu_outflow, bcu_inflow, bcu_noslip]
        return bdc
    
    def Dirichlet0Boundaries():
        """
        Defines boundaries
        """
        # Define boundary condition
        bcu_noslip = df.DirichletBC(W.sub(0), df.Constant((0, 0)), walls_boundary)
        bcu_inflow = df.DirichletBC(W.sub(0), df.Constant((0, 0)), inflow_boundary)  # parabolic input
        bcu_outflow = df.DirichletBC(W.sub(0), df.Constant((0, 0)), outflow_boundary)  # parabolic input

        bdc = [bcu_outflow, bcu_inflow, bcu_noslip]
        return bdc

    bc = DirichletBoundaries()
    bc0 = Dirichlet0Boundaries()

    up = df.TrialFunction(W)
    u,p = df.split(up)
    vq = df.TestFunction(W)
    v,q = df.split(vq)

    #Define Constant Functions and utility
    kappamin = df.Constant(1)
    kappamax = df.Constant(10)
    penalty = df.Constant(5)
    #Approximate alpha
    alpha = df.Function(K)
    #Define kappa
    kappa = kappamax + (kappamin-kappamax)*alpha*(1+penalty)/(alpha+penalty)

    #Define problem
    F = (
            -df.inner(p,df.div(v))*df.dx + #Pressure Gradient term
            df.inner(df.grad(u),df.grad(v))*df.dx + #Diffusion Term
            kappa*df.inner(u,v)*df.dx + #topology term
            df.inner(q,df.div(u))*df.dx # Divergence-Free Condition Term
        )

    up_ = df.Function(W)#Output function
    F = df.action(F,up_)#apply it to F

    # ### pyMOR wrapping
    from pymor.bindings.fenics import FenicsVectorSpace, FenicsOperator, FenicsVisualizer
    from pymor.models.basic import StationaryModel
    from pymor.operators.constructions import VectorOperator

    def alpha_setter(mu):
                alpha.vector().set_local(mu['alpha'])

    space = FenicsVectorSpace(W)
    op = FenicsOperator(F, space, space, up_,dirichlet_bcs=bc,
                        parameter_setter=alpha_setter,
                        parameters={'alpha':len(alpha.vector()[:])},
                        solver_options={'inverse': {'type': 'newton', 'rtol': 1e-6}})

    rhs = VectorOperator(op.range.zeros())

    fom = StationaryModel(op, rhs,
                          visualizer=FenicsVisualizer(space))
    return fom,K,W,bc,bc0


if __name__ == '__main__':
    main(2,16,True)

matrix and matrix_2 are not equal.
What error did I make ?

Edit :
I updated the title and adapted the code to 'handle' the boundary conditions.
Removing the boundary conditions from both the rom and the fom does make it so that matrix and matrix_2 are equal. However I am unsure as to how the boundary conditions are to be implemented here.
Using

for bc in bcs:
      bc.apply(up.vector())

(and similar for vq) should work theoretically, but due to orthonormalization, does not.

[Frostyant]

I solved the issue as follows.

First you have to add in the missing contributions from the boundaries manually via

def add_diagonal_contributions(A_ij,basis_i,basis_j,bcs):
    for bc in bcs:
        dofs = bc.get_boundary_values().keys()
        for dof in dofs:
            A_ij += basis_i[dof]*basis_j[dof]
    return A_ij

Then you have to apply 0-boundary conditions to v and compute the integral as above via assemble.
This process simulates setting a row of the matrix to [0,...,0,1,0,...,0] and the corresponding right hand-side term to the boundary value when you consider the discretized matrix equation Au=b.

In my tests it was slower than just doing it directly via pymor at higher resolutions (though faster at low resolutions and if you have a lot of basis functions).