FDM: heatmap.html con due grafici indipendenti, statico e striscia animata

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

CLAUDE.md

@@ -6,16 +6,22 @@ This file provides guidance to Claude Code (claude.ai/code) when working with co
 
 Write all git commit messages in Italian.
 
+## Scope of Work
+
+**Modifica solo il modulo `fdm/`.** Ignora qualsiasi richiesta riguardante PINN (`model.py`, `engine.py`, `visualizer.py`, `app.py` radice). Se una richiesta coinvolge file PINN, avvisa l'utente e non apportare modifiche.
+
 ## Project Overview
 
-**Heat Equation PINN** — A Physics-Informed Neural Network that solves the 1D time-varying heat equation with physical boundary conditions:
+**Heat Equation PINN** — A Physics-Informed Neural Network that solves the 1D time-varying heat equation with an internal point heat source:
 
 ```
-∂T/∂t = α ∂²T/∂x²    x ∈ [0, L],  t ∈ [0, T_END]
+∂T/∂t = α ∂²T/∂x²  +  (α/k) Q(t) δ(x − X_SRC)    x ∈ [0, L],  t ∈ [0, T_END]
 ```
 
+where `Q(t) = Q_VAL` if `t ≥ T_STEP` else `0`.
+
 - **IC:** `T(x, 0) = T0` (uniform)
-- **BC x=0:** Neumann — heat flux step: `−k ∂T/∂x = Q(t)` where `Q = Q_VAL` if `t ≥ T_STEP` else `0`
+- **BC x=0:** Robin — convection: `−k ∂T/∂x = h (T − T_AMB)`
 - **BC x=L:** Robin — convection: `−k ∂T/∂x = h (T − T_AMB)`
 
 No experimental data is needed. A `fdm/` module provides a reference numerical solution (FTCS explicit scheme) used for evaluation and visualization comparison.
@@ -59,7 +65,7 @@ model.py           ← HeatPINN + heat_pinn_loss()
 engine.py          ← data sampling, Adam+L-BFGS training, evaluation vs FDM, visualization call
 visualizer.py      ← PINN vs FDM: heatmap, animated T(x), time-series at fixed points
 fdm/
-  solver.py        ← FTCS explicit scheme, ghost-cell Neumann, explicit Robin
+  solver.py        ← FTCS explicit scheme, Robin on both ends, point source at X_SRC
   visualizer.py    ← same 3 plot types for FDM-only output
   app.py           ← FDM CLI menu
 ```
@@ -74,11 +80,11 @@ T = T_AMB + (Q_VAL · L / K) · net(x, t)
 ```
 This keeps `net` outputs in `[0, 1]` range and ensures gradients `∂T/∂x` are O(1) for the network to learn. Do not remove this scaling.
 
-`heat_pinn_loss()` normalizes all three loss terms to O(1) using `T_char = Q_VAL·L/K` and `grad_char = (Q_VAL/K)²`. Changing physical parameters in `config.py` does not require re-tuning loss weights.
+`heat_pinn_loss()` normalizes all four loss terms to O(1) using `T_char = Q_VAL·L/K` and `grad_char = (Q_VAL/K)²`. The PDE residual includes the Gaussian-smoothed source term (σ=0.02) as a continuous approximation to δ(x − X_SRC). Changing physical parameters in `config.py` does not require re-tuning loss weights.
 
 ### Training (`engine.py`)
 
-`prepare_data()` samples collocation points with **deliberate clustering**: extra points near `x=0` (steep Neumann gradient) and around `t=T_STEP` (flux step discontinuity). Increasing `N_f` / `N_bc` here is the first lever to pull if accuracy is low.
+`prepare_data()` samples collocation points with **deliberate clustering**: extra points near `x=X_SRC` (steep gradient at source) and around `t=T_STEP` (flux step discontinuity). Increasing `N_f` / `N_bc` here is the first lever to pull if accuracy is low.
 
 `train_model()` runs **Adam first, then L-BFGS fine-tuning**. L-BFGS uses a closure that captures loss components in `_last` dict (avoids calling `heat_pinn_loss` outside an active grad context).
 
@@ -87,8 +93,8 @@ This keeps `net` outputs in `[0, 1]` range and ensures gradients `∂T/∂x` are
 ### FDM Solver (`fdm/solver.py`)
 
 Returns `(T_matrix[NX, NT], x_vals, t_vals)`. Uses:
-- Ghost cell for Neumann: `T_ghost = T[1] + 2·dx·Q(t)/k`
-- Explicit Robin at x=L: `T[N] = (T[N−1] + dx·h/k·T_amb) / (1 + dx·h/k)` — uses `T_cur[-2]`, not `T_new[-2]`
+- Robin BC on both ends: `T[0] = (T[1] + robin_coeff·T_amb) / (1 + robin_coeff)`
+- Point source injected at node `i_src = argmin|x - X_SRC|` after BCs: `T[i_src] += Q·α·dt/(k·dx)`
 - CFL check at startup (warns, does not crash)
 
 ### Loss Scaling Notes

app.py

@@ -5,7 +5,7 @@ import engine
 def print_header():
     print("=" * 55)
     print("      Heat Equation PINN  —  ∂T/∂t = α ∂²T/∂x²")
-    print("      Neumann BC (x=0) + Robin BC (x=L)")
+    print("      Robin BC (x=0, x=L) + point source @ X_SRC")
     print("=" * 55)
 
 

clear.sh

@@ -0,0 +1,72 @@
+#!/usr/bin/env bash
+set -euo pipefail
+
+RESULTS_DIR="$(dirname "$0")/results/fdm"
+
+if [[ ! -d "$RESULTS_DIR" ]]; then
+    echo "Nessuna cartella results/fdm trovata."
+    exit 0
+fi
+
+mapfile -t RUNS < <(find "$RESULTS_DIR" -mindepth 1 -maxdepth 1 -type d | sort)
+
+if [[ ${#RUNS[@]} -eq 0 ]]; then
+    echo "Nessun risultato da cancellare."
+    exit 0
+fi
+
+echo "Risultati disponibili:"
+for i in "${!RUNS[@]}"; do
+    printf "  %2d. %s\n" "$((i+1))" "$(basename "${RUNS[$i]}")"
+done
+echo ""
+echo "a. Cancella tutti"
+echo "s. Selezione manuale"
+echo "0. Annulla"
+echo ""
+read -rp "Scelta: " MODE
+
+case "$MODE" in
+    a|A)
+        read -rp "Confermi la cancellazione di ${#RUNS[@]} cartelle? [s/N] " CONFIRM
+        if [[ "${CONFIRM,,}" == "s" ]]; then
+            rm -rf "${RUNS[@]}"
+            echo "Cancellati ${#RUNS[@]} risultati."
+        else
+            echo "Annullato."
+        fi
+        ;;
+    s|S)
+        read -rp "Numeri da cancellare (es. 1 3 5): " -a CHOICES
+        TO_DELETE=()
+        for N in "${CHOICES[@]}"; do
+            if [[ "$N" =~ ^[0-9]+$ ]] && (( N >= 1 && N <= ${#RUNS[@]} )); then
+                TO_DELETE+=("${RUNS[$((N-1))]}")
+            else
+                echo "Indice non valido: $N (ignorato)"
+            fi
+        done
+        if [[ ${#TO_DELETE[@]} -eq 0 ]]; then
+            echo "Nessuna selezione valida."
+            exit 0
+        fi
+        echo "Verranno cancellati:"
+        for D in "${TO_DELETE[@]}"; do
+            echo "  - $(basename "$D")"
+        done
+        read -rp "Confermi? [s/N] " CONFIRM
+        if [[ "${CONFIRM,,}" == "s" ]]; then
+            rm -rf "${TO_DELETE[@]}"
+            echo "Cancellati ${#TO_DELETE[@]} risultati."
+        else
+            echo "Annullato."
+        fi
+        ;;
+    0|"")
+        echo "Annullato."
+        ;;
+    *)
+        echo "Scelta non valida."
+        exit 1
+        ;;
+esac

engine.py

@@ -43,16 +43,17 @@ def prepare_data(N_f=4000, N_ic=400, N_bc=400):
     x_f = torch.rand(N_f, device=device) * config.L
     t_f = torch.rand(N_f, device=device) * config.T_END
 
-    # Extra points near x=0 (steep Neumann gradient) and t=T_STEP (flux step)
+    # Extra points near X_SRC (steep gradient at source) and t=T_STEP (flux step)
     n_extra = N_f // 4
-    x_near0 = torch.rand(n_extra, device=device) * config.L * 0.1          # x in [0, 0.1]
-    t_near0 = torch.rand(n_extra, device=device) * config.T_END
+    x_near_src = config.X_SRC + (torch.rand(n_extra, device=device) - 0.5) * config.L * 0.1
+    x_near_src = x_near_src.clamp(0, config.L)
+    t_near_src = torch.rand(n_extra, device=device) * config.T_END
     x_step  = torch.rand(n_extra, device=device) * config.L
     t_step  = config.T_STEP + (torch.rand(n_extra, device=device) - 0.5) * 0.1  # t near T_STEP
     t_step  = t_step.clamp(0, config.T_END)
 
-    x_f = torch.cat([x_f, x_near0, x_step])
-    t_f = torch.cat([t_f, t_near0, t_step])
+    x_f = torch.cat([x_f, x_near_src, x_step])
+    t_f = torch.cat([t_f, t_near_src, t_step])
 
     x_ic = torch.rand(N_ic, device=device) * config.L
     t_bc = torch.rand(N_bc, device=device) * config.T_END
@@ -163,10 +164,11 @@ def evaluate_model(data):
     from fdm.solver import solve as fdm_solve
     T_fdm, _, _ = fdm_solve()
 
-    # Downsample FDM time axis (NX=100, NT=5000) to match PINN grid (100x100)
+    # Downsample FDM grid (NX, NT) to match PINN prediction grid (100x100)
     nx, nt = T_pred.shape
+    x_indices = np.linspace(0, T_fdm.shape[0] - 1, nx, dtype=int)
     t_indices = np.linspace(0, T_fdm.shape[1] - 1, nt, dtype=int)
-    T_fdm_ds = T_fdm[:, t_indices]  # (100, 100)
+    T_fdm_ds = T_fdm[np.ix_(x_indices, t_indices)]  # (100, 100)
 
     l2_rel = np.sqrt(np.mean((T_pred - T_fdm_ds) ** 2)) / np.sqrt(np.mean(T_fdm_ds ** 2))
     max_err = np.max(np.abs(T_pred - T_fdm_ds))

fdm/app.py

@@ -24,26 +24,23 @@ def main_menu():
         print("\n" + "-" * 30)
         print("         MAIN MENU")
         print("-" * 30)
-        print("1. Risolvi e salva risultati")
-        print("2. Heatmap T(x,t)")
-        print("3. Animazione T(x) nel tempo")
-        print("4. Grafico T(t) in punti fissi")
+        print("1. Risolvi")
+        print("2. Visualizza")
         print("0. Esci")
         print("-" * 30)
 
-        choice = input("Select an option (0-4): ").strip()
+        choice = input("Select an option (0-2): ").strip()
 
         if choice == "1":
             T, x_vals, t_vals = solve()
             print(f"Soluzione completata. Shape T: {T.shape}")
             print(f"T range: [{T.min():.2f}, {T.max():.2f}] °C")
 
-        elif choice in ("2", "3", "4"):
+        elif choice == "2":
             if T is None:
                 print("Eseguire prima l'opzione 1.")
             else:
                 visualize_fdm(T, x_vals, t_vals)
-                print("Grafici salvati in animations/fdm/")
 
         elif choice == "0":
             print("Uscita.")

fdm/visualizer.py

@@ -1,5 +1,6 @@
 import sys
 import os
+from datetime import datetime
 import numpy as np
 import plotly.graph_objects as go
 
@@ -7,16 +8,6 @@ sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
 import config
 
 BASE_DIR = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
-FDM_ANIM_DIR = os.path.join(BASE_DIR, 'animations', 'fdm')
-
-
-def _next_path(base_dir, prefix, ext):
-    i = 1
-    while True:
-        path = os.path.join(base_dir, f'{prefix}_{i:03d}{ext}')
-        if not os.path.exists(path):
-            return path
-        i += 1
 
 
 def visualize_fdm(T_matrix, x_vals, t_vals):
@@ -31,55 +22,127 @@ def visualize_fdm(T_matrix, x_vals, t_vals):
     t_vals : np.ndarray, shape (NT,)
         Time values.
     """
-    os.makedirs(FDM_ANIM_DIR, exist_ok=True)
+    timestamp = datetime.now().strftime('%Y%m%d_%H%M%S')
+    out_dir = os.path.join(BASE_DIR, 'results', 'fdm', timestamp)
+    os.makedirs(out_dir, exist_ok=True)
 
     # ------------------------------------------------------------------
-    # 1. Heatmap
+    # 1. heatmap.html: plot statico T(x,t) + striscia 1D animata
+    #    Due figure Plotly indipendenti nello stesso file HTML per
+    #    evitare conflitti tra animazioni e assi multipli.
     # ------------------------------------------------------------------
-    zmax = float(np.max(np.abs(T_matrix - config.T0)))
-    # Use symmetric range centred on T0 if there is any variation,
-    # otherwise fall back to the raw range.
-    z_data = T_matrix.T  # shape (NT, NX) — rows=time, cols=space
+    zmin_val = float(np.min(T_matrix))
+    zmax_val = float(np.max(T_matrix))
 
-    if zmax > 0:
-        z_center = float(np.mean(T_matrix))
-        z_abs = float(np.max(np.abs(T_matrix - z_center)))
-        zmin_sym = z_center - z_abs
-        zmax_sym = z_center + z_abs
-    else:
-        zmin_sym = float(np.min(T_matrix))
-        zmax_sym = float(np.max(T_matrix))
+    # Subsample frames (shared with animation below)
+    n_frames = len(t_vals)
+    max_frames = 200
+    step = max(1, n_frames // max_frames)
+    frame_indices = list(range(0, n_frames, step))
 
-    fig_heatmap = go.Figure(go.Heatmap(
+    # --- Figura A: heatmap 2D statica (identica all'originale) ---
+    z_data = T_matrix.T  # (NT, NX)
+    z_center = float(np.mean(T_matrix))
+    z_abs = float(np.max(np.abs(T_matrix - z_center))) or 1.0
+    fig_static = go.Figure(go.Heatmap(
         z=z_data,
         x=x_vals,
         y=t_vals,
         colorscale='RdBu_r',
-        zmin=zmin_sym,
-        zmax=zmax_sym,
+        zmin=z_center - z_abs,
+        zmax=z_center + z_abs,
         colorbar=dict(title='T [°C]'),
     ))
-    fig_heatmap.update_layout(
+    fig_static.update_layout(
         title='FDM — Heat Equation  T(x,t)',
         xaxis_title='x [m]',
         yaxis_title='t [s]',
-        height=500,
+        height=480,
+        margin=dict(t=50, b=50, l=70, r=90),
     )
 
-    heatmap_path = _next_path(FDM_ANIM_DIR, 'heatmap', '.html')
-    fig_heatmap.write_html(heatmap_path)
+    # --- Figura B: striscia 1D animata ---
+    strip_frames = [
+        go.Frame(
+            data=[go.Heatmap(
+                z=T_matrix[:, idx].reshape(1, -1),
+                x=x_vals,
+                y=[0],
+                colorscale='Jet',
+                zmin=zmin_val,
+                zmax=zmax_val,
+                showscale=True,
+                colorbar=dict(title='T [°C]', thickness=18),
+            )],
+            name=str(idx),
+            layout=go.Layout(title_text=f't = {t_vals[idx]:.3f} s'),
+        )
+        for idx in frame_indices
+    ]
+
+    fig_strip = go.Figure(
+        data=[go.Heatmap(
+            z=T_matrix[:, frame_indices[0]].reshape(1, -1),
+            x=x_vals,
+            y=[0],
+            colorscale='Jet',
+            zmin=zmin_val,
+            zmax=zmax_val,
+            showscale=True,
+            colorbar=dict(title='T [°C]', thickness=18),
+        )],
+        frames=strip_frames,
+        layout=go.Layout(
+            title=f't = {t_vals[frame_indices[0]]:.3f} s',
+            xaxis=dict(title='x [m]'),
+            yaxis=dict(showticklabels=False, showgrid=False,
+                       zeroline=False, fixedrange=True),
+            height=300,
+            margin=dict(t=80, b=130, l=70, r=110),
+            updatemenus=[dict(
+                type='buttons',
+                showactive=False,
+                y=1.22, x=0.5, xanchor='center',
+                buttons=[
+                    dict(label='▶ Play', method='animate',
+                         args=[None, dict(frame=dict(duration=40, redraw=True),
+                                          fromcurrent=True, mode='immediate')]),
+                    dict(label='⏸ Pause', method='animate',
+                         args=[[None], dict(frame=dict(duration=0, redraw=False),
+                                            mode='immediate')]),
+                ],
+            )],
+            sliders=[dict(
+                steps=[
+                    dict(method='animate',
+                         args=[[str(idx)], dict(mode='immediate',
+                                                frame=dict(duration=0, redraw=True))],
+                         label=f'{t_vals[idx]:.2f}')
+                    for idx in frame_indices
+                ],
+                transition=dict(duration=0),
+                x=0.05, y=0, len=0.9,
+                currentvalue=dict(prefix='t = ', font=dict(size=14)),
+            )],
+        ),
+    )
+
+    # Scrivi entrambe le figure in un unico file HTML
+    html_static = fig_static.to_html(full_html=False, include_plotlyjs='cdn')
+    html_strip = fig_strip.to_html(full_html=False, include_plotlyjs=False)
+    heatmap_path = os.path.join(out_dir, 'heatmap.html')
+    with open(heatmap_path, 'w', encoding='utf-8') as _f:
+        _f.write('<!DOCTYPE html><html><head><meta charset="utf-8"></head><body>\n')
+        _f.write(html_static)
+        _f.write('\n')
+        _f.write(html_strip)
+        _f.write('\n</body></html>')
     print(f"Heatmap saved        → {heatmap_path}")
 
     # ------------------------------------------------------------------
     # 2. Animated profile T(x) evolving in time
     # ------------------------------------------------------------------
     L = config.L
-    n_frames = len(t_vals)
-
-    # Subsample frames for a manageable animation (max ~200 frames)
-    max_frames = 200
-    step = max(1, n_frames // max_frames)
-    frame_indices = list(range(0, n_frames, step))
 
     y_min = float(np.min(T_matrix)) - 1.0
     y_max = float(np.max(T_matrix)) + 1.0
@@ -168,7 +231,7 @@ def visualize_fdm(T_matrix, x_vals, t_vals):
         frames=frames,
     )
 
-    anim_path = _next_path(FDM_ANIM_DIR, 'animation', '.html')
+    anim_path = os.path.join(out_dir, 'animation.html')
     fig_anim.write_html(anim_path)
     print(f"Animation saved      → {anim_path}")
 
@@ -222,6 +285,6 @@ def visualize_fdm(T_matrix, x_vals, t_vals):
         height=480,
     )
 
-    ts_path = _next_path(FDM_ANIM_DIR, 'time_series', '.html')
+    ts_path = os.path.join(out_dir, 'time_series.html')
     fig_ts.write_html(ts_path)
     print(f"Time-series saved    → {ts_path}")

model.py

@@ -26,29 +26,32 @@ def heat_pinn_loss(model, x_f, t_f, x_ic, t_bc, w_pde=1.0, w_ic=1.0, w_bc=10.0):
     T_char = config.Q_VAL * config.L / config.K      # ~50 °C — temperature scale
     grad_char = (config.Q_VAL / config.K) ** 2        # ~2500 — gradient scale²
 
-    # PDE residual: dT/dt - alpha * d2T/dx2 = 0  (normalized by T_char/t_char)
+    # PDE residual: dT/dt - alpha * d2T/dx2 - source(x,t) = 0  (normalized by T_char/t_char)
     x_f = x_f.detach().requires_grad_(True)
     t_f = t_f.detach().requires_grad_(True)
     T_f = model(torch.stack([x_f, t_f], dim=1))
     dT_dt = torch.autograd.grad(T_f.sum(), t_f, create_graph=True)[0]
     dT_dx = torch.autograd.grad(T_f.sum(), x_f, create_graph=True)[0]
     d2T_dx2 = torch.autograd.grad(dT_dx.sum(), x_f, create_graph=True)[0]
+    sigma = 0.02
+    Q_t_f = torch.where(t_f >= config.T_STEP,
+                        torch.tensor(config.Q_VAL, device=t_f.device, dtype=t_f.dtype),
+                        torch.tensor(0.0,          device=t_f.device, dtype=t_f.dtype))
+    gauss = torch.exp(-0.5 * ((x_f - config.X_SRC) / sigma) ** 2) / (sigma * (2 * torch.pi) ** 0.5)
+    source_term = (config.ALPHA / config.K) * Q_t_f * gauss
     pde_scale = (T_char / config.T_END) ** 2 + 1e-8
-    L_pde = ((dT_dt - config.ALPHA * d2T_dx2) ** 2).mean() / pde_scale
+    L_pde = ((dT_dt - config.ALPHA * d2T_dx2 - source_term) ** 2).mean() / pde_scale
 
     # IC: T(x, 0) = T0 — normalized by T_char²
     T_ic_pred = model(torch.stack([x_ic, torch.zeros_like(x_ic)], dim=1))
     T_ic_true = torch.full_like(T_ic_pred, config.T0)
     L_ic = ((T_ic_pred - T_ic_true) ** 2).mean() / (T_char ** 2 + 1e-8)
 
-    # BC x=0: Neumann — dT/dx + Q(t)/K = 0
+    # BC x=0: Robin — dT/dx + H_CONV/K * (T(0,t) - T_AMB) = 0
     x_left = torch.zeros(t_bc.shape[0], device=t_bc.device).requires_grad_(True)
     T_left = model(torch.stack([x_left, t_bc.detach()], dim=1))
     dT_dx_left = torch.autograd.grad(T_left.sum(), x_left, create_graph=True)[0]
-    Q_t = torch.where(t_bc >= config.T_STEP,
-                      torch.tensor(config.Q_VAL, device=t_bc.device, dtype=t_bc.dtype),
-                      torch.tensor(0.0,          device=t_bc.device, dtype=t_bc.dtype))
-    L_bc_left = ((dT_dx_left + Q_t / config.K) ** 2).mean() / grad_char
+    L_bc_left = ((dT_dx_left + (config.H_CONV / config.K) * (T_left.squeeze() - config.T_AMB)) ** 2).mean() / grad_char
 
     # BC x=L: Robin — dT/dx + H_CONV/K * (T(L,t) - T_AMB) = 0
     x_right = torch.full((t_bc.shape[0],), config.L, device=t_bc.device).requires_grad_(True)

visualizer.py

@@ -11,10 +11,12 @@ ANIMATIONS_DIR = os.path.join(BASE_DIR, 'animations')
 def visualize_heat_field(T_pred, x_vals, t_vals, T_fdm):
     os.makedirs(ANIMATIONS_DIR, exist_ok=True)
 
-    # Downsample T_fdm from shape (NX, NT) to (NX, len(t_vals))
+    # Downsample T_fdm from shape (NX_fdm, NT_fdm) to match PINN grid
+    nx_pred = len(x_vals)
     nt_pred = len(t_vals)
+    x_indices = np.linspace(0, T_fdm.shape[0] - 1, nx_pred, dtype=int)
     t_indices = np.linspace(0, T_fdm.shape[1] - 1, nt_pred, dtype=int)
-    T_fdm_ds = T_fdm[:, t_indices]  # now shape (NX, nt_pred)
+    T_fdm_ds = T_fdm[np.ix_(x_indices, t_indices)]
 
     # --- Static heatmap: PINN vs FDM ---
     fig_map = make_subplots(