PSA Project

class: titlepage numbers

.title[Advanced Scientific Programming - Project :notebook:]

.subtitle[`PSA` - N. Dubray - ENSIIE - 2018]

.wauto[
```C++
#include <Python.h>
#include "solver.h"

static PyObject * calc_system(PyObject *self, PyObject *args)
{
  const char *string;
  int val;

  if (!PyArg_ParseTuple(args, "s", &string)) return NULL;

  Solver s(string);
  val = s.getVal();

  return Py_BuildValue("i", val);
}
```
]

.footnote[
[:book:](../index.html)
]

---

layout: true
class: animated fadeIn middle numbers

.footnote[
`PSA` - N. Dubray - ENSIIE - 2018 - [:book:](../index.html)
]

---

# Courses / project

## Schedule
.hcenter[
| Time          | 04/13         | 04/27         | 05/04                  | 05/18         | 05/25         | 06/01                  |
|:-------------:|:-------------:|:-------------:|:----------------------:|:-------------:|:-------------:|:----------------------:|
| 9:00  - 10:45 | **Course #0** | **Course #2** | **Course #4**          | **Course #5** | **Course #7** | **Course #9**          |
| 11:00 - 12:45 | Project       | Project       | Project                | Project       | Project       | Project                |
| 14:00 - 15:45 | **Course #1** | **Course #3** | .redb[Presentation #0] | **Course #6** | **Course #8** | .redb[Presentation #1] |
| 16:00 - 17:45 | Project       | Project       | .redb[Presentation #0] | Project       | Project       | .redb[Presentation #1] |
]

.vspace[]

.noflex[
.rightf.w30[
.hcenter[![](images/slides.png)]
.hcenter[[source: xkcd](https://xkcd.com/365)]
]

## Presentations
* slides [~20min] + discussion / questions [~10min]
* .redb[presentation #0: weight 1.0] (mid-project)
* .redb[presentation #1: weight 2.0] (final)

.vspace[]

## Teamwork
* 2 students per team
* **mandatory** `git` use
]

---

# Project rules

.column.middle.grow[
## Rules:
* an "`AUTHORS`" file **must** be present and contain your name(s)
* the source **must** be written in `C++11` and `Python2`
* the source **must** be documented with `Doxygen`
* the compilation chain **must** use `GNU Make`
* the presentation **must** be at the location "`pres/index.html`"
* using branches is OK but **only master will be checked**
* `Python2` bindings **are mandatory** and **must** be at the location "`bindings/`"
* the presentation **must** use `remark.js`
* unit tests **are mandatory** and **must pass**

## :moneybag: Bonuses:
* no warning at {compilation, bindings, tests, doc}
* project is implementing other methods, or is going further than specified
* source code is well written (comments, easy to read, etc...)
* commits are small and explicit
]

---

# Evaluation

## What will be evaluated at the mid-project presentation (.redb[weight 1.0])
* the presentation quality
* the teamwork
* the progress status of the project

## What will be evaluated at the final presentation (.redb[weight 2.0])
* the presentation quality
* the teamwork
* the **final** status of the project
* **the obtained results**
* **the use of the development tools**
* **unit tests**

---

# Project description

## Main objective
:arrow_right: **write and use a 2D-FD solver for the time-dependent non-relativistic Schrödinger equation**

## Details
* **write the solver** in a mix of `Python2` and `C++11` using `Armadillo`
* visualize the time-evolution of a wave packet with `Paraview` or a custom `VTK` program
* get some **real-time monitoring information** of a run using an API (in `Python2`)
* implement an **automatic restart mechanism** (in `Python2`)
* store the results in a **mongoDB database** (in `Python2`)
* calculate some **special physical cases** (diffraction, tunneling effect, HO solutions, etc...)
* use `git`, `doxygen`, `make`, `cxxtest`, `swig`, `remark.js`, `VTK` / `Paraview`

.hcenter.w40[
<video controls="yes" loop="yes" preload="auto" width="100%" height="auto" data-setup="{}" loop>
<source src="videos/tutu.webm" type="video/webm" />
</video>
]

---

# Project - formalism

## Time-dependent non-relativistic Schrödinger equation
Evolution equation:
`$$i\hbar\frac{\partial}{\partial t}\psi(x,y,t) = \hat{\mathcal{H}}_{(x,y)}\psi(x,y,t)$$`
with the 2D-Hamiltonian operator defined as
`$$\hat{\mathcal{H}}_{(x,y)}\equiv \frac{\hat{p}_{(x)}^2}{2m} + \frac{\hat{p}_{(y)}^2}{2m} + \hat{V}(x,y),$$`
and the momentum operators defined as
`$$\forall u \in \{x,y\},\hspace{5mm}\hat{p}_{(u)}\equiv -i\hbar\frac{\partial}{\partial u}.$$`
Constants:
* `\(\hbar\)` is the **reduced Planck constant**
`$$\hbar \equiv 6.582119514\times10^{−22}\,\textrm{MeV.s}$$`
* `\(m\)` is the **neutron mass**
`$$m\equiv 939.5654133\,\textrm{MeV/c}^2$$`
* `\(i\)` is the **imaginary unit** (`\(i^2=-1\)`)

---

# Project - fields

## Wave function `\(\psi(x,y,t)\)`
* is a complex field
* depends on time
* must be initialized to
`$$\psi_0(x,y) \equiv \psi(x,y,t=0),$$`
* must stay **normalized**
`$$\forall t, \iint \psi^*(x,y,t).\psi(x,y,t)\,dx\,dy = \iint \psi_0^*(x,y).\psi_0(x,y)\,dx\,dy = 1.$$`

## Potential `\(\hat{V}(x,y)\)`
* is a real field
* does not depend on time

---

# Project - solver

## Numerical method
* Finite Difference Method
* **must implement FTCS, BTCS and CTCS schemes**
* **scheme type specified at execution**

## Discretized domain
* regular mesh
* **number of points specified at execution**
* **physical span specified at execution**

## Time domain
* fixed time step
* **time step specified at execution**

## Boundary conditions
* the wave function vanishes outside the discretized domain

## Fields
* **initial wave function specified at execution**
* **potential field specified at execution**

---

# Project - architecture

.mermaid.fshadow.w60.hcenter[
    graph LR
    A("Field generator") --> B
    B("Solver") --> C("Monitor")
    style B fill:#090 
]

## Field generator
* must be written in **Python2**
* generate the initial wave function
* generate the potential field from an image file, a formula, etc...

## Solver
* must be written in **Python2** and **C++**
* **calculate the dynamical propagation**
* must be able to be **stopped and restarted automatically**
* must implement an **API** to check the progress of the calculation
* must store results in a **mongoDB database**

## Monitor
* must be written in **Python2**
* connect to the solver **or** to a mongoDB database
* generate plots / animations **or** export to *paraview*

---

# Project - special physical cases

## 2D-HO solutions
* check stationarity of some HO solutions
* check periodicity of combinations of HO solutions

## Young's slits
* check the interference pattern created by a wave packet through two slits in a potential wall

## Quantum tunneling
* check the partial crossing of a barrier higher in energy than the cinetic energy of a wave packet

## Other cases
* be creative...

---

# Project - 2D-HO solutions

## Resulting 2D-HO Schrödinger equation
`$$\hat{\mathcal{H}}_{(x,y)}^{\textrm{HO}}\psi(x,y) = E\psi(x,y)$$`
with the 2D-HO Hamiltonian operator defined as
`$$\hat{\mathcal{H}}_{(x,y)}^{\textrm{HO}}\equiv \frac{\hat{p}_{(x)}^2}{2m} + \frac{\hat{p}_{(y)}^2}{2m} + \frac{1}{2}m\omega^2(\hat{x}^2+\hat{y}^2),$$`
and the momentum operators defined as
`$$\forall u \in \{x,y\},\hspace{5mm}\hat{p}_{(u)}\equiv -i\hbar\frac{\partial}{\partial u}.$$`

## Solutions
**One can show** that solutions take the form
`$$\forall n_x\ge 0,\, \forall n_y\ge 0,\, \psi_{n_x,n_y}(x,y)=\psi_{n_x}(x).\psi_{n_y}(y),$$`
with
`$$\forall u \in \{x,y\},\hspace{5mm}\psi_{n_u}(u) = \frac{1}{\sqrt{2^n_u n_u!}}\left(\frac{m\omega}{\pi\hbar}\right)^{1/4}e^{-\frac{m\omega u^2}{2\hbar}}H_{n_u}\left(\sqrt{\frac{m\omega}{\hbar}} . u\right).$$`
The corresponding eigenvalues are given by
`$$E_{n_x, n_y}=\hbar \omega \left(n_x+n_y + 1\right).$$`

---

# Project - Young's slits

## History
* "double-slit" experiment by **Thomas Young** in 1801
* demonstrated the validity of the **wave theory of light**
* can also show the concept of **wave-particle duality** when performed with several particles

.hcenter.w60[![](images/double_slit.svg)]
.hcenter[[source: wikipedia](https://en.wikipedia.org/wiki/Double-slit_experiment#/media/File:Double-slit.svg)]

## Result
* an interference pattern appears on the screen

---

# Project - quantum tunneling

.noflex[
.rightf.w30[![](images/square_barrier.png).hcenter[[source: wikipedia](https://en.wikipedia.org/wiki/Rectangular_potential_barrier#/media/File:Square_potential.png)]]

## Setup
* a particle with **constant** energy `\(E \gt 0\)`
* a rectangular potential barrier of height `\(V_0 \gt 0\)`

## Classical mechanics
* the outcome depends only on the barrier's height, **not its shape**
* if `\(E \gt V_0\)`, the particule **always** crosses, **never** bounces
* if `\(E \lt V_0\)`, the particule **never** crosses, **always** bounces

## Quantum mechanics
* the outcome **depends on the barrier's height and shape**
* for whatever `\(E\)` and `\(V_0\)` values, the particule **may** cross, **may** bounce
* even if `\(E \gt V_0\)`, the particule **may** bounce
* even if `\(E \lt V_0\)`, the particule **may** cross :arrow_right: **quantum tunneling effect**
]

---

# Project - results presentation

:arrow_right: How to present a calculation result ?

## Input
* show the space and time discretization
* show the numerical scheme
* show the initial wave packet
* show the potential

## Output
* show the evolution with time of the total norm of the wave packet
* show the evolution with time of the wave packet

## To go further
* show projections and/or slices of the wave packet if needed
* show the (un)stability of the results with respect to the time step, the mesh size, etc...

.w70.alert.hcenter[:warning: a plot without a label and a unit for each axis is **worthless** :warning:]

---

# Project - Derivation of the FTCS scheme

## Explicit method (FTCS)

* forward difference at time `\(t_n\)`
* second-order central difference at position `\(x_i\)`

## Evolution equation
`$$i\hbar\frac{\partial}{\partial t}\psi(x,y,t) = \hat{\mathcal{H}}_{(x,y)}\psi(x,y,t)$$`
with the 2D-Hamiltonian operator defined as
`$$\hat{\mathcal{H}}_{(x,y)}\equiv \frac{\hat{p}_{(x)}^2}{2m} + \frac{\hat{p}_{(y)}^2}{2m} + \hat{V}(x,y),$$`
and the momentum operators defined as
`$$\forall u \in \{x,y\},\hspace{5mm}\hat{p}_{(u)}\equiv -i\hbar\frac{\partial}{\partial u}.$$`

---

# Project - Derivation of the BTCS scheme

## Implicit method (BTCS)

* backward difference at time `\(t_{n+1}\)`
* second-order central difference at position `\(x_i\)`

## Evolution equation
`$$i\hbar\frac{\partial}{\partial t}\psi(x,y,t) = \hat{\mathcal{H}}_{(x,y)}\psi(x,y,t)$$`
with the 2D-Hamiltonian operator defined as
`$$\hat{\mathcal{H}}_{(x,y)}\equiv \frac{\hat{p}_{(x)}^2}{2m} + \frac{\hat{p}_{(y)}^2}{2m} + \hat{V}(x,y),$$`
and the momentum operators defined as
`$$\forall u \in \{x,y\},\hspace{5mm}\hat{p}_{(u)}\equiv -i\hbar\frac{\partial}{\partial u}.$$`

---

# Project - Derivation of the CTCS scheme

## Crank-Nicolson method (CTCS)

* central difference at time `\(t_{n+\frac{1}{2}}\)`
* second-order central difference at position `\(x_i\)`

## Evolution equation
`$$i\hbar\frac{\partial}{\partial t}\psi(x,y,t) = \hat{\mathcal{H}}_{(x,y)}\psi(x,y,t)$$`
with the 2D-Hamiltonian operator defined as
`$$\hat{\mathcal{H}}_{(x,y)}\equiv \frac{\hat{p}_{(x)}^2}{2m} + \frac{\hat{p}_{(y)}^2}{2m} + \hat{V}(x,y),$$`
and the momentum operators defined as
`$$\forall u \in \{x,y\},\hspace{5mm}\hat{p}_{(u)}\equiv -i\hbar\frac{\partial}{\partial u}.$$`

---

# Project - Derivation of the Crank-Nicolson scheme

## Solution (to be demonstrated)

`$$\begin{eqnarray}\left(\frac{2i\hbar}{dt}-V(x,y)-\frac{\hbar^2}{m.dx^2}-\frac{\hbar^2}{m.dy^2}\right)\psi^{(t+dt)} &=&
\frac{-\hbar^2}{2m.dx^2}\left(\psi_{x+dx}+\psi_{x-dx}\right)\\
&-& \frac{\hbar^2}{2m.dy^2}\left(\psi_{y+dy}+\psi_{y-dy}\right)\\
&-& \frac{\hbar^2}{2m.dx^2}\left(\psi^{(t+dt)}_{x+dx}+\psi^{(t+dt)}_{x-dx}\right)\\
&-& \frac{\hbar^2}{2m.dy^2}\left(\psi^{(t+dt)}_{y+dy}+\psi^{(t+dt)}_{y-dy}\right)\\
&+& \left(V(x,y) + \frac{2i\hbar}{dt}+\frac{\hbar^2}{m.dx^2}+\frac{\hbar^2}{m.dy^2}\right)\psi
\end{eqnarray}$$`