<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 09 - Advection Dominated problem
Keywords
Step1:  3. Affine decomposition
For this problem the affine decomposition is straightforward
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 (Lagrange P2)
Step4:  4.3. Allocate an object of the AdvectionDominated 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 dolfin import *
from rbnics import *
from problems import *
from reduction_methods import *
Explanation: Tutorial 09 - Advection Dominated problem
Keywords: POD-Galerkin method, SUPG
1. Introduction
This tutorial addresses the POD-Galerkin method to the advection dominated worked problem in a two-dimensional domain $\Omega=(0,1)^2$ shown below:
<img src="data/advection_dominated_1.png" />
We introduce a stabilization technique such as $\textit{Streamline/Upwind Petrov-Galerkin}$ (SUPG) able to reduce the numerical oscillations on the approximation of the solution of parametrized advection-diffusion problem:
$$
-\varepsilon(\boldsymbol{\mu})\,\Delta u(\boldsymbol{\mu})+\beta(\boldsymbol{\mu})\cdot\nabla u(\boldsymbol{\mu})=f(\boldsymbol{\mu})\quad\text{on }\,\Omega(\boldsymbol{\mu}),
$$
where $\beta(\boldsymbol{\mu})$ and $\varepsilon(\boldsymbol{\mu})$ represent the advection and the diffusion term, respectively.
For this problem, we consider on parameter $\mu$, thus $P=1$. It is related to the Péclet number:
$$
\mathbb{P}e_K(\boldsymbol{\mu})(x):=\frac{|\beta(\boldsymbol{\mu})(x)| h_K}{2\,\varepsilon(\boldsymbol{\mu})(x)}\quad\forall x\in K\quad\forall\boldsymbol{\mu}\in\mathbb{P}.
$$
Here $h_K$ represents the diameter of $K\in\mathcal{T}_h$, where $\mathcal{T}_h$ indicates a triangulation of our domain $\Omega(\boldsymbol{\mu})$.
The parameter domain is thus given by 
$$
\mathbb{P}=[0, 6].
$$
In this problem we consider two approaches:
1. Offline-Online stabilized,
2. Offline-only stabilized,
in which while in the first one we apply the SUPG method both in the Offline and Online phases, in the second one only in the Offline phase it is applied.
In order to obtain a faster approximation of the problem, we pursue a model reduction by means of a POD-Galerkin reduced order method.
2. Parametrized formulation
Let $u(\boldsymbol{\mu})$ be the solution in the domain $\Omega$.
The PDE formulation of the parametrized problem is given by:
<center>for a given parameter $\mu=\boldsymbol{\mu}\in\mathbb{P}$, find $u(\boldsymbol{\mu})$ such that</center>
$$
\begin{cases}
    -\frac{1}{10\,^{\boldsymbol{\mu}}}\Delta\,u(\boldsymbol{\mu})+(1,1)\cdot\nabla u(\boldsymbol{\mu})=0 & \text{in }\Omega,\
    u(\boldsymbol{\mu}) = 0 & \text{on } \Gamma_1\cup\Gamma_2, \ 
    u(\boldsymbol{\mu}) = 1 & \text{on } \Gamma_3\cup\Gamma_4.
\end{cases}
$$
<br>
The corresponding weak formulation reads:
<center>for a given parameter $\boldsymbol{\mu}\in\mathbb{P}$, find $u(\boldsymbol{\mu})\in\mathbb{V}$ such that</center>
$$a\left(u(\boldsymbol{\mu}),v;\boldsymbol{\mu}\right)=f(v;\boldsymbol{\mu})\quad \forall v\in\mathbb{V},$$
where
the function space $\mathbb{V}$ is defined as
$$
\mathbb{V} = \left{ v \in H^1(\Omega): v|{\Gamma_1\cup\Gamma_2} = 0, v|{\Gamma_3\cup\Gamma_4} = 1\right},
$$
the parametrized bilinear form $a(\cdot, \cdot; \boldsymbol{\mu}): \mathbb{V} \times \mathbb{V} \to \mathbb{R}$ is defined by
$$a(u,v;\boldsymbol{\mu}) = \int_{\Omega} \frac{1}{10\,^{\boldsymbol{\mu}}}\nabla u \cdot \nabla v +\left(\partial_xu+\partial_yu\right)v\ d\boldsymbol{x},$$
the parametrized linear form $f(\cdot; \boldsymbol{\mu}): \mathbb{V} \to \mathbb{R}$ is defined by
$$f(v; \boldsymbol{\mu}) = \int_{\Omega} v\ d\boldsymbol{x}.$$
For the $\textit{Offline-Online stabilized}$ approach we use a different bilinear form $a_{stab}$ instead of $a$;
while in the $\textit{Offline-only stabilized}$ approach we use the the bilinear form $a_{stab}$ during the Offline phase, performing the Online Galerkin projection with respect to the bilinear form $a$,
the parametrized bilinear stabilized form $a_{stab}(\cdot, \cdot; \boldsymbol{\mu}): \mathbb{V} \times \mathbb{V} \to \mathbb{R}$ is defined by
$$a_{stab}(u,v,\boldsymbol{\mu}) = a(u,v,\boldsymbol{\mu}) + s(u,v,\boldsymbol{\mu}),$$
where
$$
\begin{align}
    a(u,v;\boldsymbol{\mu}) &= \int_{\Omega} \frac{1}{10\,^{\boldsymbol{\mu}}}\nabla u \cdot \nabla v +\left[(1,1)\cdot\nabla u\right]v\ d\boldsymbol{x},\
    s(u,v;\boldsymbol{\mu}) &= \sum_{K\in\mathcal{T}_h}\delta_K\int_K 
    \left(-\frac{1}{10\,^{\boldsymbol{\mu}}}\Delta u+(1,1)\cdot\nabla u\right)\left(\frac{h_K}{\sqrt{2}}(1,1)\cdot\nabla v\right)\ d\boldsymbol{x},
\end{align}
$$
and
* the parametrized linear form $f_{stab}(\cdot; \boldsymbol{\mu}): \mathbb{V} \to \mathbb{R}$ is defined by
$$
f_{stab}(v;\boldsymbol{\mu}) = f(v;\boldsymbol{\mu}) + r(v;\boldsymbol{\mu})
$$
where
$$
\begin{align}
    f(v;\boldsymbol{\mu}) &= \int_{\Omega} v\ d\boldsymbol{x}, \
    r(v;\boldsymbol{\mu}) &= \sum_{K\in\mathcal{T}_h}\delta_K\int_K \left(\frac{h_K}{\sqrt{2}}(1,1)\cdot\nabla v\right)\ d\boldsymbol{x}.
\end{align}
$$
End of explanation
@OnlineStabilization()
class AdvectionDominated(EllipticCoerciveProblem):
    # Default initialization of members
    def __init__(self, V, **kwargs):
        # Call the standard initialization
        EllipticCoerciveProblem.__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"]
        self.u = TrialFunction(V)
        self.v = TestFunction(V)
        self.dx = Measure("dx")(subdomain_data=subdomains)
        self.ds = Measure("ds")(subdomain_data=boundaries)
        # Store advection and forcing expressions
        self.beta = Constant((1.0, 1.0))
        self.f = Constant(1.0)
        # Store terms related to stabilization
        self.delta = 0.5
        self.h = CellDiameter(V.mesh())
    # Return custom problem name
    def name(self):
        return "AdvectionDominated1POD"
    # Return theta multiplicative terms of the affine expansion of the problem.
    def compute_theta(self, term):
        mu = self.mu
        if term == "a":
            theta_a0 = 10.0**(- mu[0])
            theta_a1 = 1.0
            if self.stabilized:
                delta = self.delta
                theta_a2 = - delta * 10.0**(- mu[0])
                theta_a3 = delta
            else:
                theta_a2 = 0.0
                theta_a3 = 0.0
            return (theta_a0, theta_a1, theta_a2, theta_a3)
        elif term == "f":
            theta_f0 = 1.0
            if self.stabilized:
                delta = self.delta
                theta_f1 = delta
            else:
                theta_f1 = 0.0
            return (theta_f0, theta_f1)
        else:
            raise ValueError("Invalid term for compute_theta().")
    # Return forms resulting from the discretization of the affine expansion of the problem operators.
    def assemble_operator(self, term):
        v = self.v
        dx = self.dx
        if term == "a":
            u = self.u
            beta = self.beta
            h = self.h
            a0 = inner(grad(u), grad(v)) * dx
            a1 = inner(beta, grad(u)) * v * dx
            a2 = inner(div(grad(u)), h * inner(beta, grad(v))) * dx
            a3 = inner(inner(beta, grad(u)), h * inner(beta, grad(v))) * dx
            return (a0, a1, a2, a3)
        elif term == "f":
            f = self.f
            beta = self.beta
            h = self.h
            f0 = f * v * dx
            f1 = inner(f, h * inner(beta, grad(v))) * dx
            return (f0, f1)
        elif term == "k":
            u = self.u
            k0 = inner(grad(u), grad(v)) * dx
            return (k0,)
        elif term == "m":
            u = self.u
            m0 = inner(u, v) * dx
            return (m0,)
        elif term == "dirichlet_bc":
            bc0 = [DirichletBC(self.V, Constant(0.0), self.boundaries, 1),
                   DirichletBC(self.V, Constant(0.0), self.boundaries, 2)]
            return (bc0,)
        elif term == "inner_product":
            u = self.u
            x0 = inner(grad(u), grad(v)) * dx
            return (x0,)
        else:
            raise ValueError("Invalid term for assemble_operator().")
Explanation: 3. Affine decomposition
For this problem the affine decomposition is straightforward:
$$a(u,v;\boldsymbol{\mu})=\underbrace{\frac{1}{10\,^{\boldsymbol{\mu}}}}{\Theta^{a}_0(\boldsymbol{\mu})}\underbrace{\int{\Omega}\nabla u \cdot \nabla v \ d\boldsymbol{x}}{a_0(u,v)} \ + \  \underbrace{1}{\Theta^{a}1(\boldsymbol{\mu})}\underbrace{\int{\Omega}\left[(1,1)\cdot\nabla u\right]v \ d\boldsymbol{x}}{a_1(u,v)},$$
$$f(v; \boldsymbol{\mu}) = \underbrace{1}{\Theta^{f}0(\boldsymbol{\mu})} \underbrace{\int{\Omega}v \ d\boldsymbol{x}}_{f_0(v)}.$$
Adding the following forms, we obtaing the affine decomposition for the stabilized approach:
$$s(u,v;\boldsymbol{\mu}) = \sum_{K\in\mathcal{T}h}\underbrace{\frac{\delta_K}{10\,^{\boldsymbol{\mu}}}}{\Theta^{s}0(\boldsymbol{\mu})}\underbrace{\int_K 
    \Delta u\left(\frac{h_K}{\sqrt{2}}(1,1)\cdot\nabla v\right)\ d\boldsymbol{x}}{s_0(u,v)} \ + \
    \sum_{K\in\mathcal{T}h}\underbrace{\delta_K}{\Theta^{s}1(\boldsymbol{\mu})}\underbrace{\int_K 
    \left((1,1)\cdot\nabla u\right)\left(\frac{h_K}{\sqrt{2}}(1,1)\cdot\nabla v\right)\ d\boldsymbol{x}}{s_1(u,v)},$$
$$r(v; \boldsymbol{\mu}) = \sum_{K\in\mathcal{T}h}\underbrace{\delta_K}{\Theta^{r}0(\boldsymbol{\mu})} \underbrace{\int_K\left(\frac{h_K}{\sqrt{2}}(1,1)\cdot\nabla v\right)\ d\boldsymbol{x}}{r_0(v)}.$$
We will implement the numerical discretization of the problem in the class
class AdvectionDominated(EllipticCoerciveProblem):
by specifying the coefficients $\Theta^{a}(\boldsymbol{\mu})$ and $\Theta^{f}_(\boldsymbol{\mu})$ in the method
def compute_theta(self, term):
and the bilinear forms $a(u, v)$ and linear forms $f_(v)$ in
def assemble_operator(self, term):
End of explanation
mesh = Mesh("data/square.xml")
subdomains = MeshFunction("size_t", mesh, "data/square_physical_region.xml")
boundaries = MeshFunction("size_t", mesh, "data/square_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
V = FunctionSpace(mesh, "Lagrange", 2)
Explanation: 4.2. Create Finite Element space (Lagrange P2)
End of explanation
problem = AdvectionDominated(V, subdomains=subdomains, boundaries=boundaries)
mu_range = [(0.0, 6.0)]
problem.set_mu_range(mu_range)
Explanation: 4.3. Allocate an object of the AdvectionDominated class
End of explanation
reduction_method = PODGalerkin(problem)
reduction_method.set_Nmax(15)
Explanation: 4.4. Prepare reduction with a POD-Galerkin method
End of explanation
reduction_method.initialize_training_set(100)
reduced_problem = reduction_method.offline()
Explanation: 4.5. Perform the offline phase
End of explanation
online_mu = (6.0, )
reduced_problem.set_mu(online_mu)
reduced_problem.solve(online_stabilization=True)
reduced_problem.export_solution(filename="online_solution_with_stabilization")
reduced_problem.export_error(filename="online_error_with_stabilization")
reduced_problem.solve(online_stabilization=False)
reduced_problem.export_solution(filename="online_solution_without_stabilization")
reduced_problem.export_error(filename="online_error_without_stabilization")
Explanation: 4.6. Perform an online solve
End of explanation
reduction_method.initialize_testing_set(100)
reduction_method.error_analysis(online_stabilization=True, filename="error_analysis_with_stabilization")
reduction_method.error_analysis(online_stabilization=False, filename="error_analysis_without_stabilization")
Explanation: 4.7. Perform an error analysis
End of explanation
reduction_method.speedup_analysis(online_stabilization=True, filename="speedup_analysis_with_stabilization")
reduction_method.speedup_analysis(online_stabilization=False, filename="speedup_analysis_without_stabilization")
Explanation: 4.8. Perform a speedup analysis
End of explanation
<END_TASK>