Question on Fenics Navier-Stokes

[Frostyant]

I am working on a slightly modified Navier-Stokes equation, and based my code on the pymor fenics_nonlinear demo. Since this is for a topology optimization of a pipe the parameter is a local density function in the Function space K, unlike the original piece of code. This should be the only major difference.

The attached piece of code generates a full order model for the navier stokes constraint at parameter location alpha_0 then solves it for the full order solution. It then uses this solution to build a reduced model with 1 basis and solves the constraint at alpha_0 again.

I am confused however why it appears as though the reduced model is inaccurate at the parameter location in which it is formed ? I expect that if I calculate the full order solution at a parameter location alpha and use this to form a reduced basis that the reduced order model be able to recover this reduced basis, or at least very accurately approximate it.

I tried switching from POD to Gram Schmidt (Gram Schmidt should just orthonomalize the basis, but at just 1 basis function this shouldn't make a difference) and using alternative reductors (I am planning on using the CoerciveRBreductor for an error estimate), but I cannot recover the full order solution.

Help explaining/debugging this would be appreciated.

from __future__ import print_function
from numpy import *
from fenics import *
from pymor.algorithms.newton import newton
from pymor.reductors.basic import StationaryRBReductor
from pymor.operators.ei import EmpiricalInterpolatedOperator
from pymor.algorithms.ei import ei_greedy
from pymor.algorithms.pod import pod
from pymor.bindings.fenics import *
from pymor.models.basic import *

###init/setup###
def init_bdc(W,K):
        """ 
        Defines boundaries
        """
        inflow  = 'near(x[0], 0) && (x[1]>3.5/5 && x[1]<4.5/5)'
        outflow = 'near(x[1], 0) && (x[0]>3.5/5 && x[0]<4.5/5)'
        walls   = 'near(x[0], 1) || near(x[1], 1) || (near(x[1],0) && (x[0]<=3.5/5 || x[0]>=4.5/5)) || (near(x[0],0) && (x[1]<=3.5/5 || x[1]>=4.5/5))'
        # Define boundary condition
        bcu_noslip  = DirichletBC(W.sub(0), Constant((0, 0)), walls)
        bcu_inflow  = DirichletBC(W.sub(0), Expression(("(x[1]-3.5/5)*(x[1]-4.5/5)*(-4*5)","0"),degree=2), inflow)#parabolar input
        bcu_outflow  = DirichletBC(W.sub(0).sub(0),Constant(0), outflow)#tangential = 0
        bcp_outflow = DirichletBC(W.sub(1), Constant(0), outflow)#0 pressu
        return [bcp_outflow,bcu_outflow,bcu_inflow,bcu_noslip]

def init_mesh(resolution):
        """
        Initializes Mesh
        """
        return UnitSquareMesh(resolution,resolution)


def init_VSpace(mesh,VectorSpaces,alpha0):
        V = VectorElement(VectorSpaces[0][0],mesh.ufl_cell(),VectorSpaces[0][1])
        Q = FiniteElement(VectorSpaces[1][0],mesh.ufl_cell(),VectorSpaces[1][1])
        K = FunctionSpace(mesh,"DG",0)#function space of piecewise constants
        Welem = MixedElement([V,Q])
        W = FunctionSpace(mesh,Welem)
        alpha = interpolate(alpha0,K) #Set or Reset alpha as a function in K
        return W,K,alpha

def init_fom(W,K,bdc,alpha,Re=0,kappamin=90,kappamax=90*10**4,penalty=1):
            """
            Sets up the fom object
            """
            #defines the fucntion which allows us to translate pymor parameter to fenics alpha function
            def alpha_setter(mu):
                alpha.vector().set_local(mu['alpha'])

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

            #Define Constant Functions and utility
            Re = Constant(Re)
            kappamin = Constant(kappamin)
            kappamax = Constant(kappamax)
            penalty = Constant(penalty)
            #Define kappa
            kappa = kappamax + (kappamin-kappamax)*alpha*(1+penalty)/(alpha+penalty)

            #Define problem
            F = (
                Re*inner(dot(grad(u),u),v)*dx + #Nonlinear Advection Term
                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
            )

            up_ = Function(W)#Output function
            F = action(F,up_)#apply it to F
            up = up_
            #End of Navier Stokes section

            #Wrap the fenics stuff into pymor
            space = FenicsVectorSpace(W)
            alpha_dict = {'alpha':len(alpha.vector()[:])}
            op = FenicsOperator(F,space,space,up,bdc,parameter_setter = alpha_setter,parameters= alpha_dict,solver_options={'inverse': {'type': 'newton', 'rtol': 1e-6}})
            rhs = VectorOperator(op.range.zeros())

            #Add objective functional as output
            #jop = self.init_objective_fom(space,up,alpha_setter,alpha_dict)

            #Setup full order model
            fom = StationaryModel(op,rhs, visualizer=FenicsVisualizer(space))
            return fom


def setup_rom(fom,Us,data_residuals,rtol=1e-7):
            """
            sets up fom based on input Us and data_resoduals
            """
            UU = fom.solution_space.empty()
            residuals = fom.solution_space.empty()
            for i, U in enumerate(Us):
                UU.append(U)#save fom solutions
                data_residual = data_residuals[i]
                residuals.append(data_residual['residuals'])#save fom residual

            dofs, cb, _ = ei_greedy(residuals, rtol=1e-7)
            ei_op = EmpiricalInterpolatedOperator(fom.operator, collateral_basis=cb, interpolation_dofs=dofs, triangular=True)
            rb, svals = pod(U, rtol=1e-7)
            fom_ei = fom.with_(operator=ei_op)
            reductor = StationaryRBReductor(fom_ei, rb)
            rom = reductor.reduce()
            rom = rom.with_(operator=rom.operator.with_(solver_options=fom.operator.solver_options))
            return rom
###END init/setup END###

### solvers ###
def solve_rom(rom,alpha):
            """
            Solve rom at current alpha location
            """
            mu = rom.parameters.parse(list(alpha.vector()[:]))
            #return self.rom.solve(mu)
            return newton(rom.operator, rom.rhs.as_vector(), rtol=1e-6,mu = mu, return_residuals=True)


def solve_fom(fom,alpha):
            """
            Uses a pymor fom model to resolve the constraint equation
            """
            mu = fom.parameters.parse(list(alpha.vector()[:]))#translate mu as pymor parameter
            return newton(fom.operator,fom.rhs.as_vector(), rtol=1e-6,mu = mu, return_residuals=True)#solve
            
###END solvers END###

### conversion functions ###

def pymor_np(U1):
        """
        from a pymor vector get an np vector
        """
        return transpose(U1.to_numpy()).reshape((len(transpose(U1.to_numpy())),))

def np_fenics(fenics_obj,np_):
        """
        applies numpy array to function
        """
        fenics_obj.vector().set_local(np_)
        return fenics_obj

def pymor_fenics(W,U1):
        """
        converts a pymor vector to a fenics vector
        """
        up = Function(W)
        up_np = pymor_np(U1) #convert pymor to np
        return np_fenics(up,up_np) #convert np to fenics

###END conversion functions END###


def main(resolution = 32,VectorSpaces=[["Lagrange",2],["Lagrange",1]],alpha0 = Constant(0)):
    """
    main demo script
    """
    mesh = init_mesh(resolution)

    W,K,alpha =  init_VSpace(mesh,VectorSpaces,alpha0)
    
    bdc =  init_bdc(W,K)

    fom = init_fom(W,K,bdc,alpha)

    U1,data1 = solve_fom(fom,alpha)

    rom = setup_rom(fom,[U1],[data1])

    U2,data2 = solve_rom(rom,alpha)

    up_fom = pymor_fenics(W,U1)

    up_rom = pymor_fenics(W,U2)
    
    print(up_fom)
    print('fom norm')
    print(assemble(inner(up_fom,up_fom)*dx))
    print('rom norm')
    print(assemble(inner(up_rom,up_rom)*dx))
    print('fom residual')
    print(data1['residual_norms'])
    print('rom residual')
    print(data2['residual_norms'])

#run main
main()

[sdrave]

Hi @Frostyant, I don't have much time to look into this right now. But a first thing you could do to understand the issue is to leave out the empirical interpolation and just use

reductor = StationaryRBReductor(fom, rb)

to check whether the issue is with the projection of the model or the operator interpolation.

[Frostyant]

Thank you for the tip, unfortunately I already tried this and it doesn't work. The result doesn't change whether I use empirical interpolation or not, which is expected given that the problem is linear as by default Re = 0. It seems as though the issue is with the projection itself, which I am not sure how to debug/fix.

The best I can think of is that the rom somehow doesn't interpret the parameter input correctly since that is the only difference with the demo/tutorial codes ?

[Frostyant]

The solution to this was actually simple. I forgot to project the reduced order model solution back to the full solution space before converting it into a fenics vector.

To fix it I simply needed setup_rom to also return reductor and then add the line U2 = reductor.reconstruct(U2) after U2,data2 = solve_rom(rom,alpha).

This gets the norms (technically the square of the norms) to agree.