# *FEniCS NB:* Laser Pulse Crossing a Crystal
### *File:*    fn_simp_06c.py — propagate the laser pulse, now with a static crystal, _i.e._ $\Delta \neq 0, \partial_t\Delta=0$

The Frantz-Nodvik equations
$$
\begin{align}
  \partial_t n + c\partial x &= \sigma c n \Delta, \\
  \partial_t \Delta          &= -\gamma \sigma c n \Delta, \\
\end{align}
$$
are a pair of coupled, first-order PDEs that describe
the propagation of a one-dimensional laser pulse across a
slab of laser gain material composed of two-level atoms.
Here, $x$ and $t$ denote respectively distance along the beam axis and
the time; $n(x,t)$ denotes photon number density in the medium;
and $\Delta(x,t)$ denotes the “population inversion”, $N_2 - N_1$,
giving the difference between the number density of atoms
in the excited state as compared to those in the ground state.
In addition, $c$ denotes the speed of light in the medium,
$\sigma$ the resonance absorption cross-section,
and $\gamma$ a factor related to the relative degeneracy
of the ground and excited states.

This problem has both spatial and temporal boundary conditions (BCs).
For the spatial BC, we apply Dirichlet to the left-hand edge,
and leave the right-hand edge free.
For the temporal BC, we specify the initial population inversion
across the domain of interest, and the initial profile of the
incident laser pulse.


```python
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np

# import os
# import sys
# from IPython.display import Image
# from IPython.core.interactiveshell import InteractiveShell
# InteractiveShell.ast_node_interactivity = "all"
```


```python
from fenics import *
from mshr import *
```

Describe the simulation: parameters, domain, boundary conditions, initial conditions.


```python
# set simulation parameters: domain, crystal, laser, time step
# -- simulation domain
Ld     = 1.00          # length of simulation domain (we force, see below, domain to [0,1])
mesh_d = 226           # mesh density
ds     = Ld / mesh_d   # size of individual mesh cells
# -- laser pulse
p_w  = (50 - 4) * ds   # laser pulse width
p_hw = p_w / 2         # laser pulse half-width
p_x0 = p_hw + 2 * ds   # laser pulse starting point
p_n0 = 1.00e02         # incident photon density
cl = 1.00              # speed of light in vacuum, β_0
cm = 163 * ds          # speed of light in crystal, β_c
# -- crystal
c_w   = 70 * ds        # crystal width
c_hw  = c_w / 2        # crystal half-width
c_x0  = 85 * ds        # crystal center
c_D0  = 1.00e06        # initial population inversion density
sigma = 2.70e-6        # resonance absorption cross-section
# -- time step    
T       = 183 * ds     # total simulation time
n_steps =  800         # number of time steps
dt      = T / n_steps  # size of time step
nip     = 100          # number of intervals between plots

# create domain: a 1D mesh on the unit interval [0,1]
mesh = UnitIntervalMesh(mesh_d)
V = FunctionSpace(mesh, "P", 2)

# define boundary functions: left and right ends
def on_L(x, on_boundary):
    return (on_boundary and near(x[0], 0.))

def on_R(x, on_boundary):
    return (on_boundary and near (x[0], 1.))
    
# specify boundary conditions
# (leaving RH boundary free)
bc_L = DirichletBC(V, Constant(0), on_L)
bc_R = DirichletBC(V, Constant(0), on_R)
bc = [bc_L, bc_R]
# bc = bc_L

# specify initial conditions
# -- define profile of incident laser pulse: half wave near left edge of domain
nph_init = Expression('abs(x[0] - p_x0) < p_hw ? p_n0 * cos(pi * (x[0] - p_x0) / p_w) : 0.',
                       degree=2, p_w=p_w, p_hw=p_hw, p_x0=p_x0, p_n0=p_n0)
# -- define profile of crystal's initial population inversion
Dcr_init = Expression('abs(x[0] - c_x0) < c_hw ? c_D0 * cos(pi * (x[0] - c_x0) / c_w) : 0.',
                    degree=2, c_w=c_w, c_hw=c_hw, c_x0=c_x0, c_D0=c_D0)
# -- initialize nph_j = nph_init
nph_j = interpolate(nph_init, V)
# -- initialize Dc_j = nph_init
Dcr_j = interpolate(Dcr_init, V)
```

Set up the variational problem.


```python
# define variational problem
nph = TrialFunction(V)
# Dcr = TrialFunction(V)
v   = TestFunction(V)
F = ( nph - nph_j + cm * dt * nph.dx(0) - sigma * cm * dt * Dcr_init * nph) * v * dx  # ∂n/∂t + c ∂n/∂x = σ c Δ n
# F = ( (Dcr - Dcr_j) / dt ) * v * dx  # ∂Δ/∂t = 0
a, L = lhs(F), rhs(F)
```

Execute simulation.


```python
# initializations
nph = Function(V)
# Dcr= Function(V)
t = 0
plt.figure(figsize=(13,8))
plt.title("population inversion density")

# time steps
for ii in range(n_steps):

    # update time
    t += dt

    # compute solution
    solve(a == L, nph, bc)

    # update plot at nip-step intervals
    if ii % nip == 0:
        plot(nph)

    # update previous solution
    nph_j.assign(nph)

# final envelope of laser pulse (if necessary)
if ii % nip != 0:
    plot(nph)

# finish plot
# plt.xlim([0.00, 0.20])
plt.ylim([0.00, 175.0])
# plt.axvline(c_x0 - c_hw,   c='b', lw=1.2)
# plt.axvline(c_x0 + c_hw,   c='b', lw=1.2)
plt.axvline(p_x0,          c='g', lw=1.2)
plt.axvline(p_x0 + cm * T, c='r', lw=1.2)
plt.axhline(p_n0 * 0.970,  c='b', lw=0.5)

plot(Dcr_j / (c_D0 / p_n0), ls='-.', c='b', lw=1.2)
plt.grid(True)
plt.show()
plt.close()
```

The following approach is useful for 2D and 3D systems,
when overlaying plots may not make sense.


```python
# define time-evolution function
def evolve():

    # initialization
    nph = Function(V)
    t = 0

    # time steps
    for ii in range(n_steps):

        # update current time
        t += dt

        # compute solution
        solve(a == L, nph, bc)

        # report solution at nip-step intervals
        if ii % nip == 0:
            yield nph

        # update previous solution
        nph_j.assign(nph)

    # final envelope of laser pulse (if necessary)
    if ii % nip == 0:
        yield nph
```


```python
n_rows = 3
n_cols = 5
fig_wd = 15
# default sizing here yields unit aspect ratio
plt.figure(figsize = (fig_wd, fig_wd * n_rows // n_cols))
plt.subplots_adjust(wspace=0.3, hspace=0.3)

idx = 0
for u in evolve():
    idx += 1
    plt.subplot(n_rows, n_cols, idx)
    plot(u)
#     plot(u, vmin=u_min, vmax=u_max)
```


```python

```

386f7479dd5c9e2de854ab9414afbcd80d3fd026
ipynb
Jupyter Notebook
examples/notebooks/amplifier/fn_simp_06c.ipynb
radiasoft/rslaser
096ebf41f56a1f34cea91a597f804d427bea34bc
Apache-2.0
examples/notebooks/amplifier/fn_simp_06c.ipynb
radiasoft/rslaser
096ebf41f56a1f34cea91a597f804d427bea34bc
Apache-2.0
2020-06-29T17:04:49.000Z
2022-03-28T21:39:28.000Z
examples/notebooks/amplifier/fn_simp_06c.ipynb
radiasoft/rslaser
096ebf41f56a1f34cea91a597f804d427bea34bc
Apache-2.0
2021-02-25T13:28:30.000Z
2021-02-25T13:28:30.000Z
Qwen/Qwen-72B
1. YES
2. YES