From 6a55427b60c6125a7b1bdbd2a848a542d26bcccf Mon Sep 17 00:00:00 2001 From: Davide Grilli Date: Mon, 3 Aug 2026 10:08:16 +0200 Subject: [PATCH] Valuta il flusso nodale sfruttando la separabilita della gaussiana L'impronta della sorgente e separabile in (x, s) e il fattore circonferenziale e costante durante il run: viene precalcolato una volta in prepara_stato_fem, mentre il fattore assiale e valutato sui soli n_x + 1 nodi distinti in x ed espanso con un prodotto esterno. Il risultato e identico (differenza max 1e-104) e il costo del termine noto scende da 108 a 23 us per passo, cioe -31% sul tempo totale di integrazione. Co-Authored-By: Claude Opus 5 (1M context) --- fem.py | 47 ++++++++++++++++++++++++++--------------------- 1 file changed, 26 insertions(+), 21 deletions(-) diff --git a/fem.py b/fem.py index f1301bf..6fefad1 100644 --- a/fem.py +++ b/fem.py @@ -107,6 +107,7 @@ def _massa_anello(indici: np.ndarray, lunghezza_segmento_m: float, n_nodi: int): def prepara_stato_fem( fascetta: dict, aria: dict, + sorgente: dict, sensore: dict, mesh_dati: dict | None = None, dt_s: float | None = None, @@ -114,7 +115,8 @@ def prepara_stato_fem( """Assembla le matrici FEM e fattorizza il sistema implicito. Restituisce lo stato costante per l'integrazione temporale: matrici, - solutore LU, coordinate nodali e indice del nodo osservato dal sensore. + solutore LU, coordinate nodali, peso circonferenziale della sorgente e + indice del nodo osservato dal sensore. """ if mesh_dati is None: mesh_dati = genera_mesh(fascetta) @@ -169,9 +171,14 @@ def prepara_stato_fem( carico_convezione = (h_esterno + h_interno) * area_nodale * T_ambiente carico_convezione += bordo_nodale * T_ambiente - # Coordinate nodali nel piano parametrico (x, s) della superficie media. - x_nodi_m = np.repeat(mesh_dati["x_nodi_m"], n_theta) - s_nodi_m = np.tile(raggio_m * mesh_dati["theta_nodi_rad"], n_x + 1) + # Peso circonferenziale dell'impronta della sorgente: dipende solo + # dall'offset y del percorso e dalla sigma dello spot, entrambi costanti + # durante il run, quindi è valutato una volta sola qui. La distanza è + # quella dell'immagine più vicina, perché la superficie è chiusa. + sigma_m = sorgente["sigma_punto_m"] + ds_sorgente = raggio_m * mesh_dati["theta_nodi_rad"] - sorgente["offset_y_percorso_m"] + ds_sorgente -= circonferenza_m * np.round(ds_sorgente / circonferenza_m) + peso_circonferenziale = np.exp(-0.5 * (ds_sorgente ** 2) / (sigma_m * sigma_m)) # Nodo osservato dal sensore: x più vicino a x_mm sul piano theta = 0. x_sensore_m = sensore["x_mm"] / 1000.0 @@ -185,8 +192,9 @@ def prepara_stato_fem( "C_su_dt": C_su_dt, "carico_convezione": carico_convezione, "solutore": splu(sp.csc_matrix(A)), - "x_nodi_m": x_nodi_m, - "s_nodi_m": s_nodi_m, + "x_nodi_m": mesh_dati["x_nodi_m"], + "peso_circonferenziale": peso_circonferenziale, + "forma_griglia": (n_x + 1, n_theta), "circonferenza_m": circonferenza_m, "area_nodale_m2": area_nodale, "T_ambiente_C": T_ambiente, @@ -205,8 +213,13 @@ def flusso_nodale_W_m2( Restituisce (x_riferimento_m, q_nodale) dove x_riferimento_m è la posizione della sorgente di indice 0. Ogni sorgente attiva contribuisce - con un'impronta gaussiana isotropa in (x, s); la distanza circonferenziale - è valutata sull'immagine più vicina, perché la superficie è chiusa. + con un'impronta gaussiana isotropa in (x, s). + + L'impronta gaussiana isotropa è separabile, q(x, s) = q_x(x) · q_s(s), e + il fattore circonferenziale q_s è costante nel tempo (precalcolato in + `prepara_stato_fem`). Qui si valuta quindi solo il fattore assiale sui + n_x + 1 nodi distinti in x e si espande con un prodotto esterno, invece di + valutare l'esponenziale su tutti i nodi della mesh. """ x_sensore_m = stato["x_sensore_m"] x_rif_iniziale = _x_riferimento_iniziale_m(sorgente, x_sensore_m) @@ -220,12 +233,7 @@ def flusso_nodale_W_m2( sigma = sorgente["sigma_punto_m"] q_picco = sorgente["flusso_termico_picco_W_m2"] * sorgente["efficienza_riscaldamento"] - circonferenza = stato["circonferenza_m"] - ds = stato["s_nodi_m"] - sorgente["offset_y_percorso_m"] - ds -= circonferenza * np.round(ds / circonferenza) - peso_circonferenziale = np.exp(-0.5 * (ds * ds) / (sigma * sigma)) - - q_nodale = np.zeros(stato["n_nodi"]) + peso_assiale = np.zeros(len(stato["x_nodi_m"])) for i in range(numero_sorgenti): x_i = x_riferimento + i * distanza @@ -236,13 +244,10 @@ def flusso_nodale_W_m2( continue dx = stato["x_nodi_m"] - x_i - q_nodale += ( - q_picco - * peso_circonferenziale - * np.exp(-0.5 * (dx * dx) / (sigma * sigma)) - ) + peso_assiale += np.exp(-0.5 * (dx * dx) / (sigma * sigma)) - return x_riferimento, q_nodale + q_nodale = q_picco * np.outer(peso_assiale, stato["peso_circonferenziale"]) + return x_riferimento, q_nodale.ravel() def passo_implicito_fem(stato: dict, T: np.ndarray, q_nodale: np.ndarray) -> np.ndarray: @@ -271,7 +276,7 @@ def simula_campo_fem( sorgente = cfg_run["sorgente"] sensore = cfg_run["sensore"] - stato = prepara_stato_fem(fascetta, aria, sensore, mesh_dati=mesh_dati) + stato = prepara_stato_fem(fascetta, aria, sorgente, sensore, mesh_dati=mesh_dati) dt = stato["dt_s"] indice_sensore = stato["indice_sensore"]