simulatore-induttori/plot_animazione.py

# Animazione della sezione della fascetta durante il passaggio delle sorgenti.
#
# Riproduce la fisica di run_0001 (stesso seed di simulate.py) e mostra:
# - il profilo di flusso termico q(x) sul lato esterno e le sorgenti in moto;
# - il campo di temperatura T(x, z) nella sezione lunghezza × spessore;
# - il sensore infrarosso e la temperatura nel punto osservato.

import random
from pathlib import Path

import matplotlib
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.animation import FuncAnimation, PillowWriter

from config import SIMULAZIONE
from materials import MATERIALI
from simulate import (
    calcola_skin_depth_m,
    configurazione_randomizzata,
    costruisci_solutore_implicito_2d,
    profilo_deposizione_z_1_m,
    profilo_flusso_incidente_W_m2,
)

# Istante di inizio dei fotogrammi mostrati (la simulazione parte comunque da 0).
T_INIZIO_ANIMAZIONE_S = 0.40

# Istante di fine dell'animazione.
T_FINE_ANIMAZIONE_S = 20

# Tempo simulato tra un fotogramma e il successivo.
DT_FRAME_S = 0.001

# Millisecondi tra i fotogrammi in riproduzione.
INTERVALLO_RIPRODUZIONE_MS = 1


def simula_campi(cfg_run: dict) -> dict:
    # Esegue la simulazione fino a T_FINE_ANIMAZIONE_S salvando, a ogni
    # fotogramma, campo di temperatura, profilo di flusso e stato del sensore.
    fascetta = cfg_run["fascetta"]
    aria = cfg_run["aria"]
    sorgente = cfg_run["sorgente"]
    sensore = cfg_run["sensore"]
    materiale = MATERIALI[fascetta["materiale"]]

    lunghezza = fascetta["lunghezza_mm"] / 1000.0
    spessore = fascetta["spessore_mm"] / 1000.0
    n_x = fascetta["n_nodi_x"]
    n_z = fascetta["n_nodi_z"]
    dx = lunghezza / n_x
    dz = spessore / n_z
    x_centri = (np.arange(n_x) + 0.5) * dx
    z_centri = (np.arange(n_z) + 0.5) * dz

    x_sensore = sensore["x_mm"] / 1000.0
    i_sensore = min(n_x - 1, max(0, int(x_sensore / dx)))

    dt = SIMULAZIONE["dt_interno_s"]

    if sorgente["skin_depth_fissa_m"] is None:
        skin_depth = calcola_skin_depth_m(materiale, sorgente["frequenza_hz"])
    else:
        skin_depth = float(sorgente["skin_depth_fissa_m"])

    solutore = costruisci_solutore_implicito_2d(
        n_x=n_x,
        n_z=n_z,
        dt_s=dt,
        dx_m=dx,
        dz_m=dz,
        materiale=materiale,
        h_esterno_W_m2K=aria["h_esterno_W_m2K"],
        h_interno_W_m2K=aria["h_interno_W_m2K"],
    )

    rho = materiale["densita_kg_m3"]
    cp = materiale["calore_specifico_J_kgK"]
    b_esterno = aria["h_esterno_W_m2K"] * dt / (rho * cp * dz)
    b_interno = aria["h_interno_W_m2K"] * dt / (rho * cp * dz)
    profilo_z = profilo_deposizione_z_1_m(z_centri, spessore, skin_depth)

    T = np.full((n_x, n_z), fascetta["temperatura_iniziale_C"], dtype=float)
    T_sensore = T[i_sensore, -1]
    tau_sensore = max(sensore["costante_tempo_s"], 1e-9)

    tempi, campi, flussi, x_riferimenti = [], [], [], []
    T_vere, T_lette = [], []

    prossimo_frame_t = 0.0
    t = 0.0
    while t <= T_FINE_ANIMAZIONE_S + 1e-12:
        x_rif, q_x = profilo_flusso_incidente_W_m2(sorgente, x_sensore, t, x_centri)

        rhs = T + (dt / (rho * cp)) * q_x[:, None] * profilo_z[None, :]
        rhs[:, 0] += b_esterno * aria["temperatura_ambiente_C"]
        rhs[:, -1] += b_interno * aria["temperatura_ambiente_C"]
        T = solutore.solve(rhs.ravel()).reshape(n_x, n_z)

        T_sensore += (T[i_sensore, -1] - T_sensore) * dt / tau_sensore

        if t + 1e-12 >= prossimo_frame_t:
            tempi.append(t)
            campi.append(T.copy())
            flussi.append(q_x.copy())
            x_riferimenti.append(x_rif)
            T_vere.append(T[i_sensore, -1])
            T_lette.append(T_sensore)
            prossimo_frame_t += DT_FRAME_S

        t += dt

    return {
        "tempi": np.array(tempi),
        "campi": campi,
        "flussi": flussi,
        "x_riferimenti": np.array(x_riferimenti),
        "T_vere": np.array(T_vere),
        "T_lette": np.array(T_lette),
        "x_centri_mm": x_centri * 1000.0,
        "spessore_mm": fascetta["spessore_mm"],
        "lunghezza_mm": fascetta["lunghezza_mm"],
        "x_sensore_mm": sensore["x_mm"],
        "sorgente": cfg_run["sorgente"],
    }


def main() -> None:
    rng = random.Random(SIMULAZIONE["seed"])
    cfg_run = configurazione_randomizzata(1, rng)
    dati = simula_campi(cfg_run)

    tempi = dati["tempi"]
    indice_inizio = int(np.searchsorted(tempi, T_INIZIO_ANIMAZIONE_S))
    n_frame = len(tempi) - indice_inizio

    lunghezza_mm = dati["lunghezza_mm"]
    spessore_mm = dati["spessore_mm"]
    x_vista_mm = (-10.0, lunghezza_mm + 10.0)

    q_max_MW = max(q.max() for q in dati["flussi"]) / 1e6
    T_max = max(c.max() for c in dati["campi"])

    sorgente = dati["sorgente"]
    numero_sorgenti = sorgente.get("numero_sorgenti", 1)
    distanza_mm = sorgente.get("distanza_sorgenti_m", 0.0) * 1000.0

    fig, (ax_flusso, ax_sezione, ax_storia) = plt.subplots(
        3, 1, figsize=(10, 8), height_ratios=[1.0, 1.6, 1.2],
        gridspec_kw={"hspace": 0.45},
    )
    fig.suptitle("Sezione della fascetta: sorgenti in transito e sensore")

    # Pannello 1: profilo di flusso sul lato esterno e posizioni delle sorgenti.
    linea_flusso, = ax_flusso.plot([], [], color="tab:red")
    marker_sorgenti, = ax_flusso.plot(
        [], [], "v", color="tab:red", markersize=10, clip_on=False
    )
    ax_flusso.annotate(
        "verso di marcia",
        xy=(0.28, 0.85), xytext=(0.55, 0.85), xycoords="axes fraction",
        arrowprops={"arrowstyle": "->", "color": "gray"},
        color="gray", va="center",
    )
    ax_flusso.set_xlim(*x_vista_mm)
    ax_flusso.set_ylim(0.0, q_max_MW * 1.25)
    ax_flusso.set_ylabel("q(x) [MW/m²]")
    ax_flusso.grid(True, alpha=0.3)
    ax_flusso.set_xticklabels([])

    # Pannello 2: campo di temperatura nella sezione (z verso il basso,
    # origine nel vertice in alto a sinistra come nel modello).
    immagine = ax_sezione.imshow(
        dati["campi"][indice_inizio].T,
        extent=(0.0, lunghezza_mm, spessore_mm, 0.0),
        aspect="auto",
        cmap="inferno",
        vmin=cfg_run["fascetta"]["temperatura_iniziale_C"],
        vmax=T_max,
        interpolation="bilinear",
    )
    ax_sezione.set_xlim(*x_vista_mm)
    ax_sezione.set_ylim(3.2 * spessore_mm, -0.6 * spessore_mm)
    ax_sezione.set_ylabel("z [mm]")
    ax_sezione.set_xlabel("x [mm]")
    # Sensore infrarosso sotto la parete interna (posizione schematica,
    # non in scala) con linea di vista tratteggiata.
    x_sens = dati["x_sensore_mm"]
    ax_sezione.plot([x_sens], [2.4 * spessore_mm], "^", color="tab:blue", markersize=12)
    ax_sezione.plot(
        [x_sens, x_sens], [1.1 * spessore_mm, 2.1 * spessore_mm],
        linestyle="--", color="tab:blue", linewidth=1,
    )
    ax_sezione.text(
        x_sens + 3, 2.4 * spessore_mm, "sensore IR",
        color="tab:blue", va="center",
    )
    # La colorbar è agganciata a tutti i pannelli per non restringere solo
    # quello della sezione, mantenendo allineati gli assi x.
    barra = fig.colorbar(
        immagine, ax=(ax_flusso, ax_sezione, ax_storia), pad=0.02, aspect=35
    )
    barra.set_label("T [°C]")

    # Pannello 3: temperatura nel punto osservato dal sensore.
    linea_vera, = ax_storia.plot([], [], label="T vera lato interno")
    linea_letta, = ax_storia.plot([], [], label="T sensore (con inerzia)")
    cursore = ax_storia.axvline(tempi[indice_inizio], color="gray", linewidth=0.8)
    ax_storia.set_xlim(0.0, T_FINE_ANIMAZIONE_S)
    ax_storia.set_ylim(15.0, max(dati["T_vere"].max(), dati["T_lette"].max()) * 1.08)
    ax_storia.set_xlabel("Tempo [s]")
    ax_storia.set_ylabel("T [°C]")
    ax_storia.legend(loc="upper left")
    ax_storia.grid(True, alpha=0.3)

    testo_tempo = ax_flusso.set_title(f"t = {tempi[indice_inizio]:.3f} s", loc="right")

    def aggiorna(frame: int):
        k = indice_inizio + frame
        t = tempi[k]

        linea_flusso.set_data(dati["x_centri_mm"], dati["flussi"][k] / 1e6)

        x_sorgenti_mm = (
            dati["x_riferimenti"][k] * 1000.0
            + np.arange(numero_sorgenti) * distanza_mm
        )
        visibili = (x_sorgenti_mm >= x_vista_mm[0]) & (x_sorgenti_mm <= x_vista_mm[1])
        marker_sorgenti.set_data(
            x_sorgenti_mm[visibili],
            np.full(int(visibili.sum()), q_max_MW * 1.12),
        )

        immagine.set_data(dati["campi"][k].T)

        linea_vera.set_data(tempi[: k + 1], dati["T_vere"][: k + 1])
        linea_letta.set_data(tempi[: k + 1], dati["T_lette"][: k + 1])
        cursore.set_xdata([t, t])
        testo_tempo.set_text(f"t = {t:.3f} s")

        return (
            linea_flusso, marker_sorgenti, immagine,
            linea_vera, linea_letta, cursore, testo_tempo,
        )

    animazione = FuncAnimation(
        fig, aggiorna, frames=n_frame, interval=INTERVALLO_RIPRODUZIONE_MS, blit=False
    )

    # Se il backend non è interattivo si salva una GIF invece di mostrare la finestra.
    if matplotlib.get_backend().lower() == "agg":
        percorso = Path("dataset") / "animazione_sezione.gif"
        animazione.save(
            percorso, writer=PillowWriter(fps=1000 // INTERVALLO_RIPRODUZIONE_MS)
        )
        print(f"Backend non interattivo: animazione salvata in {percorso}")
        return

    plt.show()


if __name__ == "__main__":
    main()