Présentation Projet PSA

class: titlepage

# Présentation projet PSA 

.title[
Solveur 2D-FD de l'équation dépendante du temps non relativiste de Schrödinger
]

.subtitle[
L. Melchior, T. Gennuso - ENSIIE - 2018
]

.row[
`$$\int_{-\infty}^\infty f(x)e^{-x^2} dx \simeq \sum_{i=0}^{n-1}w_if(x_i)$$`
]

---
layout: true
class: animated fadeIn middle numbers

.footnote[
Projet PSA - L. Melchior, T. Gennuso - ENSIIE 2018
]

---

# Flux de Travail

.mermaid.hcenter[
  graph LR
  A(Utilisateur) --JSON--> B
  B(generator) -- init --> D
  C(solveur) -- partial --> D
  D(MongoDB) -- partial --> E
  E(postproc) -- vtk --> F(Paraview)
  D -- partial --> C
]

---

class: toc
# Table des Matières 

1. Formalisation

2. Générateur

3. Solveur

4. Post-Traitement

5. Résultats

---

# Formalisation

Approximative solutions for the 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}.$$`

---

# Formalisation

##Méthode CTCS

Utilisation de l'approximation `\(\hspace10mm f'(z) \approx \frac{f(z + h_z) - f(z - h_z)}{2h_z} \)`
`$$i\hbar\frac{\partial}{\partial t}\psi(x,y,t) = \hat{\mathcal{H}}_{(x,y)}\psi(x,y,t)$$`
`$$\begin{eqnarray}\Leftrightarrow i\frac{\hbar}{h_t}(\psi(t+h_t) - \psi(t-h_t)) &=&
i\frac{\hbar}{2h_t}(\psi(t+h_t) - \psi + \psi - \psi(t-h_t)) 
\end{eqnarray}$$`


`$$\begin{eqnarray}\Leftrightarrow i\frac{\hbar}{h_t}(\psi(t+h_t) - \psi(t-h_t)) &=&
i\frac{\hbar}{2h_t}(FT + BT)
\end{eqnarray}$$`

Et en reprennant les formules déjà calculés, on obtient le résultats fournis dans l'énoncé.
---
class: top

# Générateur

```python
psi_bin = bson.binary.Binary(pickle.dumps(field_psi, protocol=2))
potential_bin = bson.binary.Binary(pickle.dumps(field_potential, protocol=2))

#Get the last (max in db) run_id, and add 1
try:
	current_run_id = client[dbname].data_psi.find_one(sort=[("run_id", -1)])["run_id"] + 1
	#-1 is for DESCENDING sort
except TypeError:
	current_run_id = 0

hash_p=hash_data(potential_bin)
hash_m=hash_data(psi_bin)
client[dbname]['data_psi'].insert_one({'mat':psi_bin, 'mat_bin_hash':hash_m, 'run_id':current_run_id, 'iteration':0})
client[dbname]['data_input'].insert_one({
                'field_potential':potential_bin,
                'field_potential_bin_hash':hash_p,
                'img_fact':json_file["img_factor"],
                'run_id':current_run_id,
                'nb_iter':json_file["nb_iter"],
                'dtime':json_file["dtime"],
                'mesh':mesh
                 })
client[dbname]['data_psi'].create_index([('run_id',1),('iteration',1)])
client[dbname]['data_input'].create_index([('run_id',1)])
 
```
--
* Sérialisation des matrices

--
* Définition du run_id

--
* Stockage dans mongoDB dans 2 collections

--
 * data_psi : {mat, hash, run_id, iteration}

--
 * data_input : toutes les données initiales 

--
* index sur ce qui identifie un document des collection

---

# Générateur

## Psi 

```python
kx=0.0
ky=70.0

def psi(x, y, tuple_parameters):
    return numpy.complex128(numpy.exp(-numpy.sqrt(x*x/2 + y*y/2))*numpy.exp(1j*(x*kx +y*ky)))
```

---
class: top

# Solveur — préparation

```python
def run(a):
    #a.FTCS()
    #a.BTCS(50)
    a.CTCS(50)
```

```python
current_run_id = int(input("Run id ? ")) 
data_input = client[dbname]['data_input'].find_one({'run_id':current_run_id})
try:
    bin_data=data_input['field_potential']
except TypeError:
    print("You have to choose a valid run_id")
    exit()

if hash_data(bin_data) != data_input['field_potential_bin_hash']:
    print('ERROR : no matching hash for potential')
    exit()

potential_u = pickle.loads(bin_data)
potential = numpy.array(potential_u,dtype=numpy.float64,order = 'F')
run_data = client[dbname]['data_psi'].find_one({'run_id':current_run_id}, sort = [("iteration", -1)])
try:
    bin_data=run_data['mat']
except TypeError:
    print("You have to choose a valid run_id")
    exit()

if hash_data(bin_data) != run_data['mat_bin_hash']:
    print('ERROR : not matching hash for psi')
    exit()
psi_u = pickle.loads(bin_data)
psi = numpy.array(psi_u,dtype=numpy.complex128,order = 'F')
```
--
* run_id choisit par l'utilisateur

--
* contrôle de l'intégrité des sauvegardes

--
* conversion numpy des données

--
* on récupère la dernière sauvegarde de la run_id demandée

---
class: top

# Solveur — calculs

Objet Solveur, et initialisation :
```python
a = schromod.Solveur()
a.setV(potential)
```

--
Initialisation de psi selon l'itération :
```python
if iteration != 0:
	a.setPsiInit(psi)
else:
	a.setPsiInit(normalize(psi,a.hx,a.hy))
```

--
Boucle de calcul :
```python
if iteration != 0:
    run(a)
    psi_bin = bson.binary.Binary(pickle.dumps(a.getPsi(), protocol=2))
    client[dbname]['data_psi'].replace_one({'run_id':current_run_id, 'iteration':iteration},
                        {'run_id':current_run_id, 'mat':psi_bin, 'iteration':iteration})

for i in range((iteration * img_factor), nb_iter):
    run(a)

    if i % img_factor == 0:
        print("psi norm at %d step : %f"%(i,norm_mesh(a.getPsi(),a.hx,a.hy)))
        iteration += 1
        psi_bin = bson.binary.Binary(pickle.dumps(a.getPsi(), protocol=2))
        hash_m=hash_data(psi_bin)
        client[dbname]['data_psi'].insert_one({'mat':psi_bin, 'mat_bin_hash':hash_m,
                                              'run_id':current_run_id, 'iteration':iteration})
```

--
* Redémarrage depuis une sauvegarde

--
* Sauvegarde régulière

--
 * serialisation

--
 * ajout à la collection mongoDB

---
class: top

# Solveur — à l'intérieur

```C++
int Solveur::FTCS()
{
    arma::mat ones(V.n_rows, V.n_cols, arma::fill::ones);
    std::complex<double> i(0.0,1.0);        /*The i value*/
    /*psi(t + ht)*/
    psi_next = (i*hb*ht)/(2.0*ma) *
            (
                (getPsiPlusX() - 2.0*psi + getPsiMinusX())/(hx*hx) +
                (getPsiPlusY() - 2.0*psi + getPsiMinusY())/(hy*hy)
            )
        + (ones - V*(i*ht)/hb) % psi;
    this->prepareIter();
    return 0;
}
```

```C++
int Solveur::CTCS(int iter)
{
	[...]
    for(j=0;j<iter;j++)
    {
		[...]
        this->prepareStep(psi_fat_tmp);
    }
    this->prepareIter();
    return 0;
}
```
idem pour BTCS

---
class: top

# Post-Traitement

Sélection et récupération des données
```python
current_run_id = int(input("Run id ? ")) 
data_input = client[dbname]['data_input'].find_one({'run_id':current_run_id})
potential_u = pickle.loads(data_input['field_potential'])
potential = numpy.array(potential_u,dtype=numpy.float64,order = 'C')
```

```python
imageToVTK("data/potential", pointData = {'V': numpy.reshape(potential,(nx,ny,1)).copy()})

#We don't take iteration == 0 cause it's not normalized and it disturb the plot
run_data_i = client[dbname]['data_psi'].find({'run_id':current_run_id, "iteration": {"$gt":0}}, sort = [("iteration", 1)])

#progress bar initialization
[... constants ...]
sys.stdout.write("%s[%s]" % (" " * bar_offset,  "⋅" * bar_width))
sys.stdout.flush()
sys.stdout.write("\b" * (bar_width+1)) # return to start of line, after '['

for i, data in enumerate(run_data_i):
    psi_u = pickle.loads(data['mat'])
    psi = numpy.array(psi_u,dtype=numpy.complex128,order = 'C')

    filename = 'data/anim_%07d' % (data['iteration'])
    file_str = filename + " created"

    neww = numpy.absolute(psi)
    neww = neww.reshape((nx,ny,1), order='C').copy()
	[... real ... imag ...]
    imageToVTK(filename, pointData = {'psi_norm': neww, 'real': real, 'imag': imag})
    sys.stdout.write("\r"*(len(file_str)+1) + file_str)

    if not i % true_width:
        sys.stdout.write('\r\033[%sC#' % (str(bar_offset + progress_offset)))
    progress_offset %= bar_width
    progress_offset += 1
    sys.stdout.flush()

sys.stdout.write("\n")
```

---

# Post-Traitement — affichage

```
$ ./postprocc
Run id ? 42
saving potential in potential.vti
data/anim_0000425 created [#################################⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅]
```

---

class: hcenter, middle

# Résultats — Fentes d'Young 

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

--- 
---

class: hcenter, middle

# Résultats — Réacteur 

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