# Inner problem

This notebook will use Dedalus to create a minimal working example of the solution to the inner problem:
\begin{align}
(\Gamma - \partial_{x}^2) u &= 0
\end{align}
where the penalty mask $\Gamma$ satisfies
\begin{align}
x &\to +\infty & \Gamma &\to 0\\
x &\to -\infty & \Gamma &\to 1
\end{align}
and the physical solution $u$ satisfies
\begin{align}
x &\to +\infty & \partial_x u &\to 1\\
x &\to -\infty & 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)
```

# Standard discontinuous mask

We will solve the problem numerically for a standard discontinuous mask, with finite boundary conditions at distance 10. 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 zero displacement from the reference solution $u_0 = x$.
It is important to align the numerical grid with discontinuities in the mask function


```python
def mask(x): return 1.0*(x<0)
```


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

xb0 = de.Chebyshev('x0',Nx[0],interval=(-10,0))
xb1 = de.Chebyshev('x1',Nx[1],interval=(0,10))
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'] = mask(x)

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

inner_solver = inner.build_solver()
inner_solver.solve()

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


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

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