<jupyter_start><jupyter_text>1D Reaction-Diffusion problem General formulationThis is a well known problem to FEM with Continuous Galerkin approximation. Spurious oscillations rise when certain (not rare) conditions are met. The general problem consists in:Find $u \in \mathbb{C}^2(\overline{\Omega})$ in a closed domain $\overline{\Omega}\subset \mathbb{R}^n\,(n=1,2,3)$  such that:\begin{equation}\left\{\begin{aligned}- &k \Delta u + \sigma u = f(x), \quad \forall x \in \Omega \\&\left. u \right|_{\partial \Omega} = g(x)\end{aligned}\right.\end{equation}where $\overline{\Omega} := \Omega \cup \partial \Omega$ denotes the closed domain, $\Omega$ is the open domain and $\partial \Omega$ is the domain's boundary (the closure). Note that $\Omega \cap \partial \Omega = \emptyset$. The solution is a scalar field such that $u: \overline{\Omega} \to \mathbb{R}$. Also, we have the "source term" that is a scalar field of the same form, $f: \overline{\Omega} \to \mathbb{R}$, but this function is a given data to the (direct) problem. The $g(x)$ is a prescribed function applied to the boundary, which is also a given data to the problem. The parameter $k$ (here a scalar field like $u$) sometimes is named as "diffusivity coefficient" related to the quantity $u$, in a physical context (transport of a quantity, look for Reynolds' Theorem for further clarifications). In some contexts, parameter $k$ must be a 2-rank tensor, as in anisotropic Porous Media flows. Analagously, we have the parameter $\sigma: \overline{\Omega} \to \mathbb{R}$, which is denoted (in some contexts) as the "reaction coefficient" related to the proportional rate that $u$ is changing due to generation/consumption. Parameters $k$ and $\sigma$ are given data to the (direct) problem too.P.S.: I will not discuss the ill problems formulation, neither others mathematical requirements to assert the well-posedness of the problem, uniqueness and so on. 1D simplificationIn some scenarios, the mathematical modeling approach can be simplified to the 1D case, which is written as follows:Find $u \in \mathbb{C}^2$ such that:\begin{equation}\left\{\begin{aligned}- &k u''(x) + \sigma u(x) = f(x), \quad \forall x \in ((x_1, x_2) \subset \mathbb{R}) \\&u(x_1) = u_1 \\&u(x_2) = u_2\end{aligned}\right.\end{equation}where $u_1, u_2 \in \mathbb{R}$ are the prescribed boundary condition values, $u'(x) \equiv \dfrac{d u}{d x}$ and thus $u''(x) \equiv \dfrac{d^2 u}{d x^2}$. For the sake of simplicity, I will omit the argument of $u$ function. Such aforementioned problems are known as Two-Point Boundary Value problems. Variational formulation and Galerkin approximationOne can find the derivation of the weak form, but I will just give its result below.Given the space of admissible solution:\begin{equation}\mathcal{U} (\overline{\Omega}) := \left\{ u \in H^1 (\overline{\Omega}) \left| \,u(x_1) = u_1, u(x_2) = u_2 \right. \right\}\end{equation}and the space of suitable variations as\begin{equation}\mathcal{V} (\overline{\Omega}) := \left\{ v \in H^1 (\overline{\Omega}) \left| \,u(x_1) = u(x_2) = 0 \right. \right\}\end{equation}The discretization of the above spaces is performed by Galerkin approximation with Lagrangean function space $\mathbb{P}$ of order $k$ defined over $\overline{\Omega}$. Thus we define here:\begin{equation}            \mathcal{S}_h^k(\overline{\Omega}) := \left\{ \varphi_h \in \mathcal{C}(\Omega^e): \left.\varphi\right|_{\Omega^e}             \in \mathbb{P}_k(\Omega^e), \forall \Omega^e \in \mathcal{T}_h \right\}\end{equation}the space of Lagrange polynomials subject to a domain partition $\mathcal{T}_h:=\cup \Omega^e \approx \overline{\Omega}$ and $\cap \Omega^e = \emptyset$ (non-overlapping elements). Thus\begin{equation}    \mathcal{U}_h := \mathcal{U} \cap \mathcal{S}_h^k \quad \text{and} \quad    \mathcal{V}_h := \mathcal{V} \cap \mathcal{S}_h^k\end{equation}we have the discretized spaces.Now, we can stated out our "discretized" Variational Formulation as:Find $u_h \in \mathcal{U}_h$ such that\begin{equation}a(u_h, v_h) = F(v_h), \quad \forall v_h \in \mathcal{V}_h\end{equation}where\begin{align}a(u, v) &:= \int_{\Omega} u' v' dx          + \int_{\Omega} \sigma u v dx\\F(v) &:= \int_{\Omega} f v dx\end{align}and the domain is $\overline{\Omega} \equiv [x_1, x_2]$.  A practical FEniCS exampleSolve this in FEniCS is very straightforward compared to classical FEM codes. The FEniCS framework provides a high-level interface which eases the pain at most.Just for learning purpose, let the given data be setted as:\begin{align*}&f(x) \equiv f = 1 & \\&u_1 = u_2 = 0 & \\&x_1 = 0, &x_2 = 1 \\&k = 10^{-8}, &\sigma = 1 \end{align*}With the source term as $f(x) = 1$, the exact solution can be easily obtained:\begin{equation}u(x) \approx\left\{ \begin{aligned}&u(x) = 1, &\text{if } x \in (x_1, x_2) \\&u(x) = 0, &\text{if } x = x_1 \text{ or } x = x_2 \end{aligned}\right.\end{equation}for small values of $k$, say $k = 10^{-8}$.So, how to solve the problem with FEniCS? We will construct the procedures stepwisely. * Importing all the libs we'll need.<jupyter_code>from fenics import *  # all the FEniCS namespace (not a recommended Python practice)
import matplotlib.pyplot as plt
from matplotlib import rc
import numpy as np<jupyter_output><empty_output><jupyter_text>* Defining the domain and related mesh<jupyter_code>x_left = 0.0
x_right = 1.0
numel = 15
mesh = IntervalMesh(numel, x_left, x_right)  # IntervalMesh(num_of_elements, inf_interval, sup_interval)
mesh_ref = IntervalMesh(100, x_left, x_right)<jupyter_output><empty_output><jupyter_text>* Setting the degree of the functions in the Continuous Galerkin method and the variation space. Additionaly, we define here a space to be employed in the projection of the reference analytical solution<jupyter_code>p = 1
V = FunctionSpace(mesh, "CG", p)  # "CG" stands for Continuous Galerkin, p is the degree
Vref = FunctionSpace(mesh_ref, "CG", p)<jupyter_output><empty_output><jupyter_text>* Defining a Python function which marks the boundaries<jupyter_code>def left(x, on_boundary):
    return x < 0+DOLFIN_EPS
def right(x, on_boundary):
    return x > 1-DOLFIN_EPS<jupyter_output><empty_output><jupyter_text>* Setting the prescribed boundary values<jupyter_code>u1, u2 = 0.0, 0.0
g_left = Constant(u1)
g_right = Constant(u2)<jupyter_output><empty_output><jupyter_text>* Now we modify the spaces, as expected<jupyter_code>bc_left = DirichletBC(V, g_left, left)
bc_right = DirichletBC(V, g_right, right)
dirichlet_condition = [bc_left, bc_right]<jupyter_output><empty_output><jupyter_text>* Here we define the source function over the domain<jupyter_code>f = Expression("1.0", domain=mesh, degree=p)<jupyter_output><empty_output><jupyter_text>* Now comes the good part. We set up the Trial and Test functions from the admissible space<jupyter_code>u_h = TrialFunction(V)
v_h = TestFunction(V)<jupyter_output><empty_output><jupyter_text>* The parameters of the problem<jupyter_code>k = Constant(1e-8)
sigma = Constant(1.0)<jupyter_output><empty_output><jupyter_text>* Then we write the unsymmetric form, very much like it is written in mathematical form<jupyter_code>a = k*inner(grad(v_h), grad(u_h))*dx + sigma*inner(u_h, v_h)*dx<jupyter_output><empty_output><jupyter_text>* Also we define the associated linear form<jupyter_code>L = f*v_h*dx<jupyter_output><empty_output><jupyter_text>* Now we declare the solution variable<jupyter_code>u_sol = Function(V)<jupyter_output><empty_output><jupyter_text>* Thus we set the discretized variational problem to be solved<jupyter_code>problem = LinearVariationalProblem(a, L, u_sol, dirichlet_condition)<jupyter_output><empty_output><jupyter_text>* And we go further and solve it with no difficulties!<jupyter_code>solver = LinearVariationalSolver(problem)
solver.solve()<jupyter_output><empty_output><jupyter_text>* But we will need to check if the solution is "good enough". So we compare with the available exact solution<jupyter_code>sol_exact = Expression(
    'x[0]<=0+tol || x[0]>=1 - tol ? 0 : 1', 
    degree=p+1,
    tol=DOLFIN_EPS
)
u_e = interpolate(sol_exact, Vref)
u_e = interpolate(sol_exact, Vref)<jupyter_output><empty_output><jupyter_text>* Now, lets plot our results!<jupyter_code>plot(u_sol, marker='x', label='Approx')
plot(u_e, label='Exact')
# Setting the font
plt.rc('text',usetex=True)
plt.rc('font', size=14)
# Plotting
plt.xlim(x_left, x_right)   # Limites do eixo x
plt.ylim(np.min(u_sol.vector().get_local()), 1.02*np.max(u_sol.vector().get_local()))  # Limites do eixo y
plt.grid(True, linestyle='--')  # Ativa o grid do grafico
plt.xlabel(r'$x$')  # Legenda do eixo x
plt.ylabel(r'$u(x)$')   # Legenda do eixo y
plt.legend(loc='best',borderpad=0.5)    # Ativa legenda no grafico e diz para se posicionar na melhor localizacao detectada
plt.show()  # Exibe o grafico em tela<jupyter_output><empty_output><jupyter_text>Something seems not good... (Batman?) Terrible, it isn't? This is due the high Damköhler number of the problem $(Da >> 1)$. The Damköhler number is conceptually defined as\begin{equation}Da := \frac{\text{reaction rate of a physical quantity}}{\text{diffusion rate of a physical quantity}}\end{equation}In Finite Element context, we evaluate the Damköhler number in an element-wise manner, because this is how we (normally) stabilize the method and vanish the spurious oscillations. The local (within each element, valid for 1D and 2D... I need to confirm for 3D) Damköhler number is denote as\begin{equation}Da_K \equiv \frac{\sigma h_k^2}{6k}\end{equation}where $h_K$ is the mesh parameter, related to some characteristic lenght of the element. For further details about what I stated above and for the method I will summarize below, check the paper "The Galerkin Gradient Least-Squares Method" (Franca and <NAME>, 1989). Highly recommend reading for who want to dive in FEM Stabilization methods.So, where all this stuff will be applied?  GGLS stabilization As commented above, this method was proposed by "The Galerkin Gradient Least-Squares Method" (Franca and <NAME>, 1989) aiming to solve singular diffusion problems in reaction dominated cases ($Da >> 1$). The principle is similar to GLS method, but control is gained in $H^1$-seminorm, while classical Galerkin (and GLS) controls over $L^2$-norm.With the modification, the variational problem is rewritten as follows:\begin{equation}B(u,v) = F(v)\end{equation}and\begin{align}&B(u,v) := a(u,v) + R(u,v) \\&F(v) := L(v) + S(v)\end{align}with\begin{align}&R(u,v) := \sum_{K}\tau_k \left(\nabla\left(\sigma u -k \Delta u\right),\nabla\left(\sigma v -k \Delta v\right)\right)_K\\&S(v) := \sum_{K}\tau_k \left(f, \nabla\left(\sigma v -k \Delta v\right) \right)_K\end{align}where\begin{equation}(u, v)_{K} := \int_{\Omega^K} u v dx\end{equation}denotes the $L^2$ inner product restricted to the element $K$.Note the presence of a stabilizing parameter $\tau_k$. By this parameter, we control how much we want to add to an element local contribution aiming to stabilize in the local formulation. This parameter is computed locally as follows:\begin{equation}\tau_K := \frac{h_K^2}{6 \sigma} \widetilde{\xi}(Da_K)\end{equation}where\begin{equation}\xi(Da_K) :=\left\{\begin{aligned}&1, &Da_K \geq 8, \\&0.064 Da_K + 0.49, &1 \leq Da_K < 8 \\&0, &Da_K < 1\end{aligned}\right.\end{equation}The above expression is an asymptotic approximation to $\widetilde{\xi}$. See the aforementioned paper for further details.  1D GGLSTo simplify to 1D, just remember that the gradient reduce to $\nabla (\bullet) \equiv \dfrac{d (\bullet)}{d x} = u'$. Thus,\begin{align}&R(u,v) := \sum_{K}\tau_k \left(\sigma u' -k u''', \sigma v' -k v''' \right)_K \\&S(v) := \sum_{K}\tau_k \left(f, \sigma u' -k v''' \right)_K\end{align}  Revisiting the FEniCS implementation Finally, we can stabilize with GGLS the previous Continuous Galerkin FEM approach. * First, we need to define the stabilizing parameter and its related quantities<jupyter_code>h = CellDiameter(mesh)
Da_k = (sigma*pow(h,2))/(6.0*k)
eps = conditional(ge(Da_k,8),1,conditional(ge(Da_k,1),0.064*Da_k+0.49,0))
tau = (eps*(pow(h,2)))/(6.0*sigma)<jupyter_output><empty_output><jupyter_text>* Now we add the stabilizing terms to the classical Galerkin formulation<jupyter_code>a += inner(grad(sigma*u_h - k*div(grad(u_h))), tau*grad(sigma*v_h - k*div(grad(v_h))))*dx
L += inner(grad(f), tau*grad(sigma*v_h - k*div(grad(v_h))))*dx<jupyter_output><empty_output><jupyter_text>* Redefining the problem<jupyter_code>u_sol = Function(V)
problem = LinearVariationalProblem(a, L, u_sol, dirichlet_condition)<jupyter_output><empty_output><jupyter_text>* Solving the stabilized problem<jupyter_code>solver = LinearVariationalSolver(problem)
solver.solve()<jupyter_output><empty_output><jupyter_text>* Now we plot the results to check<jupyter_code>plot(u_sol, marker='x', label='Approx')
plot(u_e, label='Exact')
# Configuracoes de fontes no grafico
plt.rc('text',usetex=True)
plt.rc('font', size=14)
# Plotando
plt.xlim(x_left, x_right)   # Limites do eixo x
plt.ylim(np.min(u_sol.vector().get_local()), 1.05*np.max(u_e.vector().get_local()))  # Limites do eixo y
plt.grid(True, linestyle='--')  # Ativa o grid do grafico
plt.xlabel(r'$x$')  # Legenda do eixo x
plt.ylabel(r'$u(x)$')   # Legenda do eixo y
plt.legend(loc='best',borderpad=0.5)    # Ativa legenda no grafico e diz para se posicionar na melhor localizacao detectada
plt.show()  # Exibe o grafico em tela<jupyter_output><empty_output><jupyter_text>Now things look right.  Another stabilization: Galerkin Least Squares (GLS)An alternative stabilization is the very well known method named Galerkin Least Squares, or GLS for short.This methods was introduced by Hughes, Franca and Hulbert (see the paper "A new finite element formulation for computational fluid dynamics: VIII. The Galerkin/least-squares method for advective–diffusive equations" for further details). The formulation to reaction-diffusion problem was based on a remark provide in the Franca and Do Carmo GGLS paper.It's very similar approach to GGLS method in the modification sense related to the classical Galerkin formulation. Now, we have the following terms to add:\begin{align}&R_{GLS}(u,v) := \sum_{K}\tau_k \left(\sigma u -k \Delta u, \sigma v -k \Delta v \right)_K \\&S_{GLS}(v) := \sum_{K}\tau_k \left(f, \sigma v -k \Delta v \right)_K\end{align}which trivially is extended to 1D as\begin{align}&R_{GLS}(u,v) := \sum_{K}\tau_k \left(\sigma u -k u'', \sigma v -k v''\right)_K \\&S_{GLS}(v) := \sum_{K}\tau_k \left(f, \sigma v -k v'' \right)_K\end{align}  1D GLS implementation in FEniCSFor instance, the stabilization parameter and its parameters can be pretty much the same as in GLS, just for testing purpose. Minor changes are necessary, but we will solve within a different variational problem just to compare the methods. * The classical Continuous Galerkin terms<jupyter_code>a_GLS = k*inner(grad(v_h), grad(u_h))*dx + sigma*inner(u_h, v_h)*dx
L_GLS = f*v_h*dx<jupyter_output><empty_output><jupyter_text>* Adding the GLS terms<jupyter_code>a_GLS += tau*inner((sigma*u_h - k*div(grad(u_h))), sigma*v_h - k*div(grad(v_h)))*dx
L_GLS += tau*inner(f, sigma*v_h - k*div(grad(v_h)))*dx<jupyter_output><empty_output><jupyter_text>* Now we solve the problem with GLS<jupyter_code>u_sol_gls = Function(V)
problem = LinearVariationalProblem(a_GLS, L_GLS, u_sol_gls, dirichlet_condition)
solver = LinearVariationalSolver(problem)
solver.solve()<jupyter_output><empty_output><jupyter_text>* Then we plot:<jupyter_code>plot(u_sol, marker='.', label='GGLS')
plot(u_sol_gls, marker='x', label='GLS')
plot(u_e, label='Exact')
# Configuracoes de fontes no grafico
plt.rc('text',usetex=True)
plt.rc('font', size=14)
# Plotando
plt.xlim(x_left, x_right)   # Limites do eixo x
plt.ylim(np.min(u_sol.vector().get_local()), 1.02*np.max(u_sol_gls.vector().get_local()))  # Limites do eixo y
plt.grid(True, linestyle='--')  # Ativa o grid do grafico
plt.xlabel(r'$x$')  # Legenda do eixo x
plt.ylabel(r'$u(x)$')   # Legenda do eixo y
plt.legend(loc='best',borderpad=0.5)    # Ativa legenda no grafico e diz para se posicionar na melhor localizacao detectada
plt.show()  # Exibe o grafico em tela<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13407407407407407, 'fraction_numerical': 0.009876543209876543, 'mean_word_length': 5.406089363384737, 'pattern_counts': {'":': 0, '<': 107, '<?xml version=': 0, '>': 109, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '25', 'dir': 'jupyter-structured-clean-dedup', 'id': '53335', 'license': 'MIT', 'n_tokens_mistral': 5723, 'n_tokens_neox': 5469, 'n_words': 1760, 'path': 'fem/fenics/reaction-diffusion1D.ipynb', 'repo_name': 'volpatto/quick-ipython-notebooks', 'size': '16204'}