# Heat Equation PINN

![Python](https://img.shields.io/badge/Python-3.10%2B-blue)
![PyTorch](https://img.shields.io/badge/PyTorch-2.0%2B-orange)
![Plotly](https://img.shields.io/badge/Plotly-5.15%2B-green)
![Tests](https://img.shields.io/badge/tests-42%20passing-brightgreen)

Soluzione dell'equazione del calore 1D con sorgente puntuale tramite **Physics-Informed Neural Network (PINN)**, validata contro un solver numerico **FDM** (Finite Difference Method).

---

## Indice

1. [Problema fisico](#1-problema-fisico)
2. [Approccio: PINN vs FDM](#2-approccio-pinn-vs-fdm)
3. [Struttura del progetto](#3-struttura-del-progetto)
4. [Installazione](#4-installazione)
5. [Utilizzo](#5-utilizzo)
6. [Architettura della rete neurale](#6-architettura-della-rete-neurale)
7. [Funzione di loss](#7-funzione-di-loss)
8. [Training](#8-training)
9. [Solver FDM di riferimento](#9-solver-fdm-di-riferimento)
10. [Visualizzazioni](#10-visualizzazioni)
11. [Parametri configurabili](#11-parametri-configurabili)
12. [Test](#12-test)

---

## 1. Problema fisico

Il progetto risolve l'equazione del calore 1D con una sorgente di calore puntuale interna che si attiva a un istante fissato:

```
∂T/∂t = α ∂²T/∂x²  +  (α/k) Q(t) δ(x − X_SRC)
```

dove:
- `T(x, t)` — temperatura [°C]
- `α` — diffusività termica [m²/s]
- `k` — conducibilità termica [W/m·K]
- `Q(t) = Q_VAL` se `t ≥ T_STEP`, altrimenti `0` — flusso di calore [W/m²]
- `δ(x − X_SRC)` — delta di Dirac nella posizione della sorgente
- Dominio: `x ∈ [0, L]`, `t ∈ [0, T_END]`

### Condizioni al contorno (Robin / convezione)

Entrambe le estremità della barra scambiano calore per convezione con l'ambiente:

```
x = 0:   −k ∂T/∂x = h (T − T_AMB)     →    ∂T/∂x − (h/k)(T − T_AMB) = 0
x = L:   −k ∂T/∂x = h (T − T_AMB)     →    ∂T/∂x + (h/k)(T − T_AMB) = 0
```

### Condizione iniziale

```
T(x, 0) = T₀    (temperatura uniforme)
```

### Parametri fisici di default

| Parametro    | Valore  | Unità    | Significato                              |
|--------------|---------|----------|------------------------------------------|
| `ALPHA`      | 0.01    | m²/s     | Diffusività termica                      |
| `K`          | 1.0     | W/m·K    | Conducibilità termica                    |
| `L`          | 1.0     | m        | Lunghezza della barra                    |
| `T0`         | 20.0    | °C       | Temperatura iniziale uniforme            |
| `X_SRC`      | 0.35    | m        | Posizione della sorgente di calore       |
| `Q_VAL`      | 150.0   | W/m²     | Intensità del flusso di calore           |
| `T_STEP`     | 0.2     | s        | Istante di attivazione della sorgente    |
| `H_CONV`     | 10.0    | W/m²·K   | Coefficiente convettivo alle estremità   |
| `T_AMB`      | 20.0    | °C       | Temperatura ambiente                     |
| `T_END`      | 10.0    | s        | Fine della simulazione                   |

---

## 2. Approccio: PINN vs FDM

### PINN (Physics-Informed Neural Network)

Una PINN è una rete neurale che impara a soddisfare un'equazione differenziale alle derivate parziali **senza dati sperimentali**. L'addestramento avviene minimizzando una funzione di loss che penalizza le violazioni della PDE, delle condizioni iniziali e delle condizioni al contorno in un insieme di **punti di collocazione** distribuiti nel dominio.

Vantaggi:
- Non richiede una griglia regolare
- Può essere addestrata su domini irregolari
- La soluzione è una funzione continua e differenziabile ovunque

In questo progetto la loss ha tre componenti:
1. **Residuo PDE** — la rete deve soddisfare l'equazione del calore
2. **Condizione iniziale** — la rete deve restituire `T₀` per `t = 0`
3. **Condizioni al contorno** — la rete deve soddisfare le Robin BC a `x=0` e `x=L`

### FDM (Finite Difference Method)

Il solver FDM (`fdm/`) implementa lo schema **FTCS** (Forward-Time Centered-Space) esplicito su una griglia regolare `NX × NT`. Non usa reti neurali: è una soluzione numerica classica che serve come **riferimento ad alta fedeltà** per validare la PINN.

La valutazione della PINN consiste nel calcolare l'errore relativo L2 tra la predizione della rete e la soluzione FDM interpolata sulla stessa griglia `100 × 100`.

---

## 3. Struttura del progetto

```
pinn/
├── config.py          ← tutti i parametri fisici e numerici (modificare qui)
├── model.py           ← HeatPINN (rete neurale) + heat_pinn_loss()
├── engine.py          ← campionamento dati, training, valutazione vs FDM
├── visualizer.py      ← grafici PINN vs FDM: heatmap, animazione, serie temporali
├── app.py             ← CLI interattiva per PINN
├── requirements.txt
├── pytest.ini
├── clear.sh           ← script di pulizia artefatti
├── fdm/
│   ├── solver.py      ← schema FTCS esplicito con Robin BC e sorgente puntuale
│   ├── visualizer.py  ← grafici FDM standalone
│   └── app.py         ← CLI interattiva per FDM
└── tests/
    ├── conftest.py
    ├── test_config.py
    ├── test_model.py
    ├── test_engine_data.py
    ├── test_fdm_solver.py
    ├── test_integration_pinn.py
    └── test_e2e.py
```

| File               | Responsabilità                                                          |
|--------------------|-------------------------------------------------------------------------|
| `config.py`        | Unica fonte di verità per tutti i parametri                             |
| `model.py`         | Definizione della rete e calcolo della loss fisica                      |
| `engine.py`        | Pipeline completa: campionamento → training → valutazione               |
| `visualizer.py`    | Plot interattivi HTML (PINN vs FDM)                                     |
| `app.py`           | Menu CLI per l'utente                                                   |
| `fdm/solver.py`    | Solver numerico FTCS per la soluzione di riferimento                    |
| `fdm/visualizer.py`| Plot interattivi HTML (FDM standalone)                                  |
| `fdm/app.py`       | Menu CLI per il solver FDM                                              |

---

## 4. Installazione

**Prerequisiti:** Python 3.10+, `pip`, `virtualenv` (o `venv`).

```bash
git clone <repository-url>
cd pinn

python -m venv .venv
source .venv/bin/activate        # Linux / macOS
# .venv\Scripts\activate         # Windows

pip install -r requirements.txt
```

Il progetto rileva automaticamente GPU/MPS/CPU all'avvio (vedi [engine.py — device detection](#12-dettagli-implementativi)).

---

## 5. Utilizzo

Attivare sempre il virtual environment prima di eseguire qualsiasi script:

```bash
source .venv/bin/activate
```

### PINN (`app.py`)

```bash
python app.py
```

```
1. Addestra nuovo modello
2. Valuta vs FDM (L2 error, max error)
3. Visualizza (genera 3 file HTML)
0. Esci
```

- **Opzione 1** — avvia l'addestramento (Adam + L-BFGS). Chiede il numero di epoche; premi Invio per usare il default (5000). Il modello migliore viene salvato in `models/best_heat_pinn_model.pth`. Per riaddestrare da zero: `rm models/best_heat_pinn_model.pth`.
- **Opzione 2** — carica il modello salvato, esegue il solver FDM, stampa l'errore relativo L2 e l'errore massimo assoluto.
- **Opzione 3** — genera tre file HTML interattivi in `results/pinn/<timestamp>/`.

### FDM (`fdm/app.py`)

```bash
python fdm/app.py
```

```
1. Risolvi (schema FTCS)
2. Visualizza (genera 3 file HTML)
0. Esci
```

- **Opzione 1** — esegue il solver e stampa shape della matrice e range di temperatura.
- **Opzione 2** — genera tre file HTML in `results/fdm/<timestamp>/`.

### Test

```bash
pytest tests/ -v                          # tutti i test (42)
pytest tests/test_model.py -v             # rete e loss
pytest tests/test_engine_data.py -v       # campionamento dati
pytest tests/test_fdm_solver.py -v        # solver FDM
pytest tests/test_integration_pinn.py -v  # integrazione PINN
pytest tests/test_e2e.py -v               # workflow completo
```

### Pulizia artefatti

```bash
./clear.sh    # menu interattivo per eliminare models/, results/ o entrambi
```

---

## 6. Architettura della rete neurale

`HeatPINN` ([model.py](model.py)) è una rete fully connected a 5 layer:

```
Input (x, t)
    │
[Linear 2→128, Tanh]
[Linear 128→128, Tanh]   ×3
[Linear 128→1]
    │
Output: T
```

La rete riceve le coordinate `(x, t)` e produce un unico scalare: la temperatura `T(x, t)`.

### Normalizzazione dell'input

Prima di entrare nella rete, le coordinate vengono normalizzate al range `[0, 1]`:

```python
x_norm = x / L        # x ∈ [0, 1]
t_norm = t / T_END    # t ∈ [0, 1]
```

Questo migliora il condizionamento numerico dell'ottimizzazione.

### Output scaling

La rete interna `net` non predice direttamente la temperatura: predice una **perturbazione adimensionale**. La temperatura fisica viene ricostruita con:

```python
T = T_AMB + (Q_VAL * L / K) * net(x_norm, t_norm)
```

dove `T_char = Q_VAL * L / K ≈ 150 °C` è la scala caratteristica di temperatura del problema.

**Perché questo scaling?**
- L'output della rete rimane nell'ordine di `[0, 1]`, rendendo il training più stabile.
- I gradienti `∂T/∂x` risultanti sono `O(1)` — la rete può imparare la struttura spaziale senza problemi di scala.
- Il termine di fondo `T_AMB` garantisce che la soluzione parta dalla condizione iniziale corretta anche con pesi random.

**Non rimuovere questo scaling**: senza di esso la rete deve imparare ordini di grandezza diversi tra le condizioni iniziali/al contorno e il gradiente interno, rendendo l'ottimizzazione molto più difficile.

---

## 7. Funzione di loss

`heat_pinn_loss()` ([model.py](model.py)) calcola quattro valori: `(total, L_pde, L_ic, L_bc)`.

```
total = W_PDE * L_pde + W_IC * L_ic + W_BC * L_bc
```

Ogni termine è **normalizzato automaticamente** da scale precompilate che dipendono solo dalle costanti in `config.py`. Cambiare `Q_VAL`, `K`, `H_CONV` o `L` ribilancia automaticamente la loss senza richiedere rituning manuale dei pesi.

### L_pde — Residuo PDE

Valutato su `N_F` punti di collocazione `(x_f, t_f)` distribuiti nel dominio:

```
residuo = dT/dt  −  α d²T/dx²  −  source(x, t)
```

Il termine sorgente usa un'approssimazione gaussiana al delta di Dirac (il delta non è differenziabile):

```
source(x, t) = (α/k) · Q(t) · G(x)

G(x) = exp(−0.5 · ((x − X_SRC) / σ)²) / (σ √(2π))    σ = GAUSS_SIGMA = 0.01 m
```

`Q(t)` è la funzione gradino: vale `Q_VAL` per `t ≥ T_STEP`, zero altrimenti.

I gradienti `dT/dt`, `dT/dx`, `d²T/dx²` sono calcolati con `torch.autograd.grad` — **autodifferenziazione esatta**, non differenze finite.

Normalizzazione: `_pde_scale = max((T_char/T_END)², src_peak²)` dove `src_peak` è il picco gaussiano.

### L_ic — Condizione iniziale

Valutato su `N_IC` punti `(x_ic, 0)`:

```
L_ic = mean( (T(x, 0) − T₀)² ) / T_char²
```

### L_bc — Condizioni al contorno Robin

Valutato su `N_BC` istanti temporali `t_bc`, applicato a entrambe le estremità:

```
x = 0:   ∂T/∂x(0, t) − (h/k)(T(0,t) − T_AMB) = 0
x = L:   ∂T/∂x(L, t) + (h/k)(T(L,t) − T_AMB) = 0
```

Normalizzazione: `_bc_scale = max(Q_VAL/K, H_CONV·T_char/K)²`

> **Nota sul segno:** la BC a sinistra (x=0) ha segno negativo davanti al termine convettivo perché il flusso uscente è orientato verso `−x`; a destra (x=L) il segno è positivo perché il flusso uscente è orientato verso `+x`.

---

## 8. Training

Il training è implementato in `train_model()` ([engine.py](engine.py)) e procede in due fasi.

### Fase 1: Adam

```
Ottimizzatore: Adam, LR = 1e-3
Scheduler:     ReduceLROnPlateau (factor=0.5, patience=150, min_lr=1e-6)
Early stopping: se la loss non migliora di > 1e-7 per 500 epoche consecutive
```

Il modello con la loss più bassa viene salvato a ogni miglioramento in `models/best_heat_pinn_model.pth`.

### Fase 2: L-BFGS (fine-tuning)

Al termine dell'Adam, viene caricato il miglior modello e affinato con L-BFGS:

```
Ottimizzatore: L-BFGS, LR=0.1, max_iter=50, history_size=50, strong Wolfe
Steps:         200
```

L-BFGS è un ottimizzatore di secondo ordine (quasi-Newton) particolarmente efficace nella fase finale del training PINN perché sfrutta la curvatura della loss per convergere a minimi più precisi di quanto Adam riesca a fare.

**Meccanismo closure:** L-BFGS richiede di poter rivalutare la loss più volte per iterazione. La funzione `closure()` cattura i componenti della loss in un dizionario `_last` per poterli stampare senza ricalcolare il grafo computazionale fuori dal contesto di `backward()`.

### Campionamento dei punti di collocazione

`prepare_data()` genera i punti di collocazione con **clustering deliberato** nelle zone fisicamente più complesse:

| Zona                        | Proporzione | Motivazione                                    |
|-----------------------------|-------------|------------------------------------------------|
| Uniforme `[0,L] × [0,T_END]`| 50%         | Copertura generale del dominio                 |
| Intorno a `X_SRC ± 5% L`    | 25%         | Gradiente ripido in prossimità della sorgente  |
| Intorno a `T_STEP ± 0.1 s`  | 25%         | Discontinuità temporale all'attivazione        |

Il clustering aumenta la densità di punti dove la fisica è più difficile da apprendere, senza aumentare il costo computazionale totale.

---

## 9. Solver FDM di riferimento

`fdm/solver.py` implementa lo schema **FTCS** (Forward-Time Centered-Space) esplicito.

### Schema di avanzamento temporale

```
T[i, n+1] = T[i, n]  +  r · (T[i+1,n] − 2·T[i,n] + T[i−1,n])

r = α·dt/dx²    (numero di Courant-Friedrichs-Lewy)
```

### Stabilità CFL

Condizione necessaria per la stabilità dello schema esplicito:

```
r = α·dt/dx² ≤ 0.5
```

Se la condizione è violata, il solver stampa un avvertimento ma non si blocca. Con i parametri di default (`NX=250`, `NT=15000`) la condizione è soddisfatta. Se si riducono `NX` o `NT`, verificare che `r ≤ 0.5`.

### Condizioni al contorno Robin

Le BC sono applicate a ogni passo temporale usando uno schema centrato:

```
T[0]  = (T[1]  + robin_coeff · T_AMB) / (1 + robin_coeff)   # x = 0
T[-1] = (T[-2] + robin_coeff · T_AMB) / (1 + robin_coeff)   # x = L

robin_coeff = dx · h / k
```

### Iniezione della sorgente puntuale

Dopo l'applicazione delle BC, la sorgente viene iniettata al nodo più vicino a `X_SRC`:

```
T[i_src] += Q(t) · α · dt / (k · dx)
```

---

## 10. Visualizzazioni

Tutti i plot sono **HTML interattivi** generati con Plotly (zoom, hover, slider, play/pause).

### PINN vs FDM (`visualizer.py`)

Generati con `python app.py → opzione 3`, salvati in `results/pinn/<timestamp>/`:

| File              | Contenuto                                                              |
|-------------------|------------------------------------------------------------------------|
| `heatmap.html`    | Heatmap 2D affiancate PINN vs FDM, stessa scala colori                 |
| `animation.html`  | Profilo T(x) animato nel tempo: PINN (blu continuo) vs FDM (rosso tratteggiato) |
| `comparison.html` | Serie temporali T(t) nei punti fissi `x=0`, `x=L/2`, `x=L`; linea verticale a `t=T_STEP` |

### FDM standalone (`fdm/visualizer.py`)

Generati con `python fdm/app.py → opzione 2`, salvati in `results/fdm/<timestamp>/`:

| File              | Contenuto                                                              |
|-------------------|------------------------------------------------------------------------|
| `heatmap.html`    | Heatmap 2D T(x,t) + striscia animata del profilo di temperatura        |
| `animation.html`  | Profilo T(x) animato nel tempo                                         |
| `time_series.html`| Serie temporali in 5 punti fissi: `0`, `0.25L`, `0.5L`, `0.75L`, `L` |

---

## 11. Parametri configurabili

Tutti i parametri si trovano in `config.py`. Modificare solo questo file per cambiare il problema o il training.

### Fisica

| Parametro    | Default | Unità  | Descrizione                              |
|--------------|---------|--------|------------------------------------------|
| `ALPHA`      | 0.01    | m²/s   | Diffusività termica                      |
| `K`          | 1.0     | W/m·K  | Conducibilità termica                    |
| `L`          | 1.0     | m      | Lunghezza della barra                    |
| `T0`         | 20.0    | °C     | Temperatura iniziale uniforme            |
| `X_SRC`      | 0.35    | m      | Posizione della sorgente di calore       |
| `Q_VAL`      | 150.0   | W/m²   | Intensità del flusso di calore           |
| `T_STEP`     | 0.2     | s      | Istante di attivazione della sorgente    |
| `H_CONV`     | 10.0    | W/m²·K | Coefficiente convettivo alle estremità   |
| `T_AMB`      | 20.0    | °C     | Temperatura ambiente                     |
| `T_END`      | 10.0    | s      | Fine della simulazione                   |

### Griglia FDM

| Parametro    | Default | Descrizione                                             |
|--------------|---------|---------------------------------------------------------|
| `NX`         | 250     | Nodi spaziali (aumentare per maggiore risoluzione)      |
| `NT`         | 15000   | Passi temporali (verificare `r = α·dt/dx² ≤ 0.5`)      |
| `GAUSS_SIGMA`| 0.01    | Larghezza del picco gaussiano nella loss PINN [m]       |

### Architettura PINN

| Parametro         | Default | Descrizione                                  |
|-------------------|---------|----------------------------------------------|
| `HIDDEN_SIZE`     | 128     | Neuroni per layer nascosto                   |
| `N_HIDDEN_LAYERS` | 4       | Numero di layer nascosti (totale: 5 layer)   |

### Campionamento

| Parametro | Default | Descrizione                                                        |
|-----------|---------|--------------------------------------------------------------------|
| `N_F`     | 6000    | Punti PDE (+ 50% clustering automatico vicino a X_SRC e T_STEP)   |
| `N_IC`    | 400     | Punti condizione iniziale                                          |
| `N_BC`    | 400     | Punti condizioni al contorno                                       |

### Training Adam

| Parametro       | Default | Descrizione                                         |
|-----------------|---------|-----------------------------------------------------|
| `EPOCHS`        | 5000    | Epoche massime                                      |
| `PATIENCE`      | 500     | Early stopping: epoche senza miglioramento          |
| `LR_ADAM`       | 1e-3    | Learning rate iniziale                              |
| `SCHED_FACTOR`  | 0.5     | Fattore di riduzione LR (ReduceLROnPlateau)         |
| `SCHED_PATIENCE`| 150     | Patience per la riduzione LR                        |
| `SCHED_MIN_LR`  | 1e-6    | Learning rate minimo                                |

### Fine-tuning L-BFGS

| Parametro     | Default | Descrizione              |
|---------------|---------|--------------------------|
| `LR_LBFGS`   | 0.1     | Learning rate L-BFGS     |
| `LBFGS_STEPS` | 200     | Numero di step L-BFGS    |

### Pesi della loss

| Parametro | Default | Descrizione                    |
|-----------|---------|--------------------------------|
| `W_PDE`   | 10.0    | Peso residuo PDE               |
| `W_IC`    | 1.0     | Peso condizione iniziale        |
| `W_BC`    | 5.0     | Peso condizioni al contorno    |

> **Se la loss diverge:** verificare che `T_char = Q_VAL * L / K` non sia vicino a zero. Questo valore è la scala caratteristica di temperatura usata per normalizzare tutti i termini.

---

## 12. Test

La test suite è in `tests/` e conta **42 test** organizzati in tre livelli:

| File                         | Tipo        | Cosa testa                                        |
|------------------------------|-------------|---------------------------------------------------|
| `test_config.py`             | Unit        | Validità e coerenza dei parametri in `config.py`  |
| `test_model.py`              | Unit        | Shape output, finitezza, loss components          |
| `test_engine_data.py`        | Unit        | Campionamento e clustering dei punti              |
| `test_fdm_solver.py`         | Unit        | Griglia, CFL, shape output del solver FDM         |
| `test_integration_pinn.py`   | Integration | Caricamento modello, griglia predizione, pipeline loss |
| `test_e2e.py`                | End-to-end  | Workflow completo: train → evaluate → visualize   |

```bash
pytest tests/ -v              # run tutti i test
pytest tests/ -k "model"      # solo test con "model" nel nome
pytest tests/ --tb=short      # traceback breve in caso di fallimento
```

---

## Dettagli implementativi

### Normalizzazione automatica della loss

Le scale sono precompilate una sola volta a import time in `model.py`:

```python
_T_char    = Q_VAL * L / K                      # ~150 °C
_src_peak  = ALPHA * Q_VAL / (K * GAUSS_SIGMA * sqrt(2π))
_pde_scale = max((_T_char / T_END)², _src_peak²) + 1e-8
_bc_scale  = max(Q_VAL / K, H_CONV * _T_char / K) ** 2
```

Dividere ogni termine per la sua scala porta tutti i contributi a `O(1)`, rendendo i pesi `W_PDE`, `W_IC`, `W_BC` interpretabili come importanza relativa piuttosto che come fattori di scala assoluti.

### Rilevamento device

`engine.py` seleziona automaticamente il device più performante disponibile:

```python
CUDA  →  MPS (Apple Silicon)  →  CPU
```

Include un test di funzionamento effettivo della GPU prima di usarla, per evitare fallimenti silenziosi su driver incompleti.

### Closure L-BFGS

L-BFGS richiede una funzione `closure()` che esegue `zero_grad`, forward pass, `backward`, e restituisce la loss. I componenti della loss vengono catturati in un dizionario `_last` per permettere il logging a ogni step senza ricalcolare il grafo fuori dal contesto `backward`.

### Subsampling delle animazioni FDM

Se `NT > 200`, il visualizer FDM campiona ogni `n`-esimo frame (`n = NT // 200`) per mantenere le animazioni HTML leggere e fluide.