From b663a89abde727910c48918bd61fb5ab8ff7b026 Mon Sep 17 00:00:00 2001 From: Davide Grilli Date: Thu, 14 May 2026 12:07:14 +0200 Subject: [PATCH] Allinea PINN alla fisica FDM: sorgente interna e BC Robin bilaterali MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - model.py: aggiunge termine sorgente Gaussiana (σ=0.02) nella PDE loss per approssimare δ(x − X_SRC); sostituisce BC Neumann a x=0 con Robin - engine.py: clustering collocation vicino X_SRC anziché x=0; downsample FDM su entrambi gli assi spaziale e temporale in evaluate_model() - visualizer.py: downsample FDM su entrambi gli assi prima del plot - app.py: aggiorna header con fisica corrente - CLAUDE.md: aggiorna PDE, BC e note architetturali Co-Authored-By: Claude Sonnet 4.6 --- CLAUDE.md | 18 ++++++++++-------- app.py | 2 +- engine.py | 16 +++++++++------- model.py | 17 ++++++++++------- visualizer.py | 6 ++++-- 5 files changed, 34 insertions(+), 25 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index cfbaab8..bc36411 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -8,14 +8,16 @@ Write all git commit messages in Italian. ## 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 +61,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 +76,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 +89,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 diff --git a/app.py b/app.py index 44bbada..62cef94 100644 --- a/app.py +++ b/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) diff --git a/engine.py b/engine.py index 1afb783..f21bdc7 100644 --- a/engine.py +++ b/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)) diff --git a/model.py b/model.py index 6600c19..75152f7 100644 --- a/model.py +++ b/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) diff --git a/visualizer.py b/visualizer.py index be21687..eb6940e 100644 --- a/visualizer.py +++ b/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(