<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Tutorial 17 - Navier Stokes equations
Keywords
Step1:  3. Affine Decomposition
Step2:  4. Main program
4.1. Read the mesh for this problem
The mesh was generated by the data/generate_mesh.ipynb notebook.
Step3:  4.2. Create Finite Element Space (Taylor-Hood P2-P1)
Step4:  4.3. Allocate an object of the NavierStokes class
Step5:  4.4. Prepare reduction with a POD-Galerkin method
Step6:  4.5. Perform the offline phase
Step7:  4.6. Perform an online solve
Step8:  4.7. Perform an error analysis
Step9:  4.8. Perform a speedup analysis
<ASSISTANT_TASK:>
Python Code:
from ufl import transpose
from dolfin import *
from rbnics import *
Explanation: Tutorial 17 - Navier Stokes equations
Keywords: exact parametrized functions, supremizer operator
1. Introduction
In this tutorial, we will study the Navier-Stokes equations over the two-dimensional backward-facing step domain $\Omega$ shown below:
<img src="data/backward_facing_step.png" width="80%"/>
A Poiseuille flow profile is imposed on the inlet boundary, and a no-flow (zero velocity) condition is imposed on the walls. A homogeneous Neumann condition of the Cauchy stress tensor is applied at the outflow boundary.
The inflow velocity boundary condition is characterized by $$\boldsymbol{u}(\boldsymbol{x};\mu)=\mu\bigg {\frac{1}{2.25}(x_1-2)(5-x_1),0\bigg } \quad \forall \boldsymbol{x}=(x_0,x_1) \in \Omega$$ 
This problem is characterized by one parameter $\mu$, which characterizes the inlet velocity. The range of $\mu$ is the following $$\mu \in [1.0, 80.0].$$ 
Thus, the parameter domain is $$\mathbb{P}=[1.0,80.0].$$
In order to be able to compare the interpolation method (DEIM) used to solve this problem, we pursue an exact model reduction by means of a POD-Galerkin reduced order method.
2. Parametrized formulation
Let $\boldsymbol{u}(\mu)$ be the velocity vector and $p(\mu)$ be the pressure in the domain $\Omega$.
We will directly provide a weak formulation for this problem: <center>for a given parameter $\mu \in \mathbb{P},$ find $u(\mu) \in \mathbb{V}(\mu), \; p \in\mathbb{M}$ such that </center>
<center>
    $
    \begin{cases}
        \nu \int_{\Omega_o} \nabla \boldsymbol{u_o} : \nabla \boldsymbol{v_o} \ d\Omega + \int_{\Omega_o} [(\boldsymbol{u_o} \cdot \nabla) \boldsymbol{u_o}] \cdot \boldsymbol{v_o} \ d\Omega - \int_{\Omega_o} p_o \nabla \cdot \boldsymbol{v_o} \ d\Omega = \int_{\Omega_o} \boldsymbol{f_o} \cdot \boldsymbol{v_o} \ d\Omega, \quad \forall \boldsymbol{v_o} \in\mathbb{V_o},  \
        \int_{\Omega_o} q_o \nabla \cdot \boldsymbol{u_o} \ d\Omega = 0, \quad \forall q_o \in\mathbb{M_o}
    \end{cases}
    $
</center>
where
$\nu$ represents kinematic viscosity
the functional space $\mathbb{V}(\mu)$ is defined as $\mathbb{V}=[H^1_{\Gamma_{wall}}(\Omega)]^2$
the functional space $\mathbb{M}(\mu)$ is defined as $\mathbb{M}=L^2(\Omega)$
Since this problem utilizes mixed finite element discretization with the velocity and pressure as solution variables, the inf-sup condition is necessary for the well posedness of this problem. Thus, the supremizer operator $T^{\mu}: \mathbb{M}_h \rightarrow \mathbb{V}_h$ will be used.
End of explanation
@ExactParametrizedFunctions()
@PullBackFormsToReferenceDomain()
@ShapeParametrization(
    ("x[0]", "x[1]"),  # subdomain 1
    ("x[0]", "mu[1] / 2. * x[1] + (2. - mu[1])"),  # subdomain 2
)
class NavierStokes(NavierStokesProblem):
    # Default initialization of members
    def __init__(self, V, **kwargs):
        # Call the standard initialization
        NavierStokesProblem.__init__(self, V, **kwargs)
        # ... and also store FEniCS data structures for assembly
        assert "subdomains" in kwargs
        assert "boundaries" in kwargs
        self.subdomains, self.boundaries = kwargs["subdomains"], kwargs["boundaries"]
        dup = TrialFunction(V)
        (self.du, self.dp) = split(dup)
        (self.u, _) = split(self._solution)
        vq = TestFunction(V)
        (self.v, self.q) = split(vq)
        self.dx = Measure("dx")(subdomain_data=self.subdomains)
        self.ds = Measure("ds")(subdomain_data=self.boundaries)
        # ... as well as forcing terms and inlet velocity
        self.inlet = Expression(("1. / 2.25 * (x[1] - 2) * (5 - x[1])", "0."), degree=2)
        self.f = Constant((0.0, 0.0))
        self.g = Constant(0.0)
        # Customize nonlinear solver parameters
        self._nonlinear_solver_parameters.update({
            "linear_solver": "mumps",
            "maximum_iterations": 20,
            "report": True
        })
    # Return custom problem name
    def name(self):
        return "NavierStokesExact2"
    # Return theta multiplicative terms of the affine expansion of the problem.
    @compute_theta_for_derivatives
    @compute_theta_for_supremizers
    def compute_theta(self, term):
        mu = self.mu
        if term == "a":
            theta_a0 = 1.
            return (theta_a0,)
        elif term in ("b", "bt"):
            theta_b0 = 1.
            return (theta_b0,)
        elif term == "c":
            theta_c0 = 1.
            return (theta_c0,)
        elif term == "f":
            theta_f0 = 1.
            return (theta_f0,)
        elif term == "g":
            theta_g0 = 1.
            return (theta_g0,)
        elif term == "dirichlet_bc_u":
            theta_bc00 = mu[0]
            return (theta_bc00,)
        else:
            raise ValueError("Invalid term for compute_theta().")
    # Return forms resulting from the discretization of the affine expansion of the problem operators.
    @assemble_operator_for_derivatives
    @assemble_operator_for_supremizers
    def assemble_operator(self, term):
        dx = self.dx
        if term == "a":
            u = self.du
            v = self.v
            a0 = inner(grad(u) + transpose(grad(u)), grad(v)) * dx
            return (a0,)
        elif term == "b":
            u = self.du
            q = self.q
            b0 = - q * div(u) * dx
            return (b0,)
        elif term == "bt":
            p = self.dp
            v = self.v
            bt0 = - p * div(v) * dx
            return (bt0,)
        elif term == "c":
            u = self.u
            v = self.v
            c0 = inner(grad(u) * u, v) * dx
            return (c0,)
        elif term == "f":
            v = self.v
            f0 = inner(self.f, v) * dx
            return (f0,)
        elif term == "g":
            q = self.q
            g0 = self.g * q * dx
            return (g0,)
        elif term == "dirichlet_bc_u":
            bc0 = [DirichletBC(self.V.sub(0), self.inlet, self.boundaries, 1),
                   DirichletBC(self.V.sub(0), Constant((0.0, 0.0)), self.boundaries, 2)]
            return (bc0,)
        elif term == "inner_product_u":
            u = self.du
            v = self.v
            x0 = inner(grad(u), grad(v)) * dx
            return (x0,)
        elif term == "inner_product_p":
            p = self.dp
            q = self.q
            x0 = inner(p, q) * dx
            return (x0,)
        else:
            raise ValueError("Invalid term for assemble_operator().")
# Customize the resulting reduced problem
@CustomizeReducedProblemFor(NavierStokesProblem)
def CustomizeReducedNavierStokes(ReducedNavierStokes_Base):
    class ReducedNavierStokes(ReducedNavierStokes_Base):
        def __init__(self, truth_problem, **kwargs):
            ReducedNavierStokes_Base.__init__(self, truth_problem, **kwargs)
            self._nonlinear_solver_parameters.update({
                "report": True,
                "line_search": "wolfe"
            })
    return ReducedNavierStokes
Explanation: 3. Affine Decomposition
End of explanation
mesh = Mesh("data/backward_facing_step.xml")
subdomains = MeshFunction("size_t", mesh, "data/backward_facing_step_physical_region.xml")
boundaries = MeshFunction("size_t", mesh, "data/backward_facing_step_facet_region.xml")
Explanation: 4. Main program
4.1. Read the mesh for this problem
The mesh was generated by the data/generate_mesh.ipynb notebook.
End of explanation
element_u = VectorElement("Lagrange", mesh.ufl_cell(), 2)
element_p = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
element = MixedElement(element_u, element_p)
V = FunctionSpace(mesh, element, components=[["u", "s"], "p"])
Explanation: 4.2. Create Finite Element Space (Taylor-Hood P2-P1)
End of explanation
problem = NavierStokes(V, subdomains=subdomains, boundaries=boundaries)
mu_range = [(1.0, 70.0), (1.5, 2.5)]
problem.set_mu_range(mu_range)
Explanation: 4.3. Allocate an object of the NavierStokes class
End of explanation
reduction_method = PODGalerkin(problem)
reduction_method.set_Nmax(20)
Explanation: 4.4. Prepare reduction with a POD-Galerkin method
End of explanation
lifting_mu = (1.0, 2.0)
problem.set_mu(lifting_mu)
reduction_method.initialize_training_set(100, sampling=EquispacedDistribution())
reduced_problem = reduction_method.offline()
Explanation: 4.5. Perform the offline phase
End of explanation
online_mu = (70.0, 1.5)
reduced_problem.set_mu(online_mu)
reduced_solution = reduced_problem.solve()
plot(reduced_solution, reduced_problem=reduced_problem, component="u")
plot(reduced_solution, reduced_problem=reduced_problem, component="p")
Explanation: 4.6. Perform an online solve
End of explanation
reduction_method.initialize_testing_set(36, sampling=EquispacedDistribution())
reduction_method.error_analysis()
Explanation: 4.7. Perform an error analysis
End of explanation
reduction_method.speedup_analysis()
Explanation: 4.8. Perform a speedup analysis
End of explanation
<END_TASK>