# Poiseuille flow

This notebook will use Dedalus to create a minimal working example of the solution to planar Poiseuille flow
\begin{align}
\partial_{x}^2 u - \frac{1}{\varepsilon^2}\Gamma u &= -2
\end{align}
where the penalty mask $\Gamma$ can be chosen by shifting and rescaling optimal masks $\Gamma^*$ for the inner problem
\begin{align}
\Gamma(x) &= \Gamma^*\left(\frac{x - \varepsilon \ell}{\varepsilon \delta}\right)
\end{align}
and the physical solution $u$ satisfies the boundary conditions
\begin{align}
x &= 1 & u &= 0\\
x &\to -\infty & \partial_x u &\to 0
\end{align}

# Imports


```python
import numpy as np
import matplotlib.pyplot as plt
import dedalus.public as de

from matplotlib import rc
rc('font',**{'family':'serif','serif':['Computer Modern Roman']})
rc('text', usetex=True)
```

# Optimized discontinuous mask

We will solve the problem numerically for a standard discontinuous mask, with finite boundary conditions at distance -5. We enforce Robin boundary conditions within the solid to match onto the analytical exponential behaviour.

One can change the mask function to calculate different solutions.
Masks chosen according to the optimal criterion will achieve $\mathcal{O}(\varepsilon^2)$ error from the reference solution $u_0 = x(1-x)$.


```python
# Calculate penalized solution
ϵ = 0.1
Nx = [128,128]

xb0 = de.Chebyshev('x0',Nx[0],interval=(-5,0))
xb1 = de.Chebyshev('x1',Nx[1],interval=(0,1))
xbasis = de.Compound('x',[xb0,xb1])
domain = de.Domain([xbasis], grid_dtype=np.float64)
x = xbasis.grid(*domain.dealias)

Γ = domain.new_field(name='Γ',scales=domain.dealias)
Γ['g'] = 1.*(x<ϵ)

poiseuille = de.LBVP(domain, variables=['u','ux'])
poiseuille.meta[:]['x']['dirichlet'] = True
poiseuille.parameters['Γ'] = Γ
poiseuille.parameters['ε'] = ϵ
poiseuille.add_equation("dx(ux) - (Γ/ε**2)*u = -2")
poiseuille.add_equation("ux - dx(u) = 0")
poiseuille.add_bc("left(ux) = 0")
poiseuille.add_bc("right(u) = 0")

poiseuille_solver = poiseuille.build_solver()
poiseuille_solver.solve()

u, ux = poiseuille_solver.state['u'], poiseuille_solver.state['ux']
```


```python
# Plot penalized and reference solutions
fig, ax = plt.subplots()
ax.plot(x,u['g'],color='C1',label='Penalized')
ax.plot(x[x>0],x[x>0]*(1-x[x>0]),'k--',label='Reference')
ax.fill_between(x[x<0],0,10,color='lightgray')
ax.set(aspect=1,xlim=[-1,1],ylim=[0,1],xlabel='$x$',ylabel='$u$')
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.legend()
```

94839aff525c6bc5f96a0b75fe87d12db97a0451
ipynb
Jupyter Notebook
poiseuille-flow.ipynb
ericwhester/volume-penalty-code
66a1745daeae2ad71bda0bc9299c8b8271a9871f
MIT
2020-03-14T19:40:40.000Z
2022-03-18T03:02:33.000Z
poiseuille-flow.ipynb
ericwhester/volume-penalty-code
66a1745daeae2ad71bda0bc9299c8b8271a9871f
MIT
poiseuille-flow.ipynb
ericwhester/volume-penalty-code
66a1745daeae2ad71bda0bc9299c8b8271a9871f
MIT
Qwen/Qwen-72B
1. YES
2. YES