commit 0f4438cc7fbf7c03645765d221ebda34339a232f Author: Davide Grilli Date: Wed May 13 21:21:53 2026 +0200 Commit iniziale diff --git a/.gitattributes b/.gitattributes new file mode 100644 index 0000000..1ff0c42 --- /dev/null +++ b/.gitattributes @@ -0,0 +1,63 @@ +############################################################################### +# Set default behavior to automatically normalize line endings. +############################################################################### +* text=auto + +############################################################################### +# Set default behavior for command prompt diff. +# +# This is need for earlier builds of msysgit that does not have it on by +# default for csharp files. +# Note: This is only used by command line +############################################################################### +#*.cs diff=csharp + +############################################################################### +# Set the merge driver for project and solution files +# +# Merging from the command prompt will add diff markers to the files if there +# are conflicts (Merging from VS is not affected by the settings below, in VS +# the diff markers are never inserted). Diff markers may cause the following +# file extensions to fail to load in VS. An alternative would be to treat +# these files as binary and thus will always conflict and require user +# intervention with every merge. To do so, just uncomment the entries below +############################################################################### +#*.sln merge=binary +#*.csproj merge=binary +#*.vbproj merge=binary +#*.vcxproj merge=binary +#*.vcproj merge=binary +#*.dbproj merge=binary +#*.fsproj merge=binary +#*.lsproj merge=binary +#*.wixproj merge=binary +#*.modelproj merge=binary +#*.sqlproj merge=binary +#*.wwaproj merge=binary + +############################################################################### +# behavior for image files +# +# image files are treated as binary by default. +############################################################################### +#*.jpg binary +#*.png binary +#*.gif binary + +############################################################################### +# diff behavior for common document formats +# +# Convert binary document formats to text before diffing them. This feature +# is only available from the command line. Turn it on by uncommenting the +# entries below. +############################################################################### +#*.doc diff=astextplain +#*.DOC diff=astextplain +#*.docx diff=astextplain +#*.DOCX diff=astextplain +#*.dot diff=astextplain +#*.DOT diff=astextplain +#*.pdf diff=astextplain +#*.PDF diff=astextplain +#*.rtf diff=astextplain +#*.RTF diff=astextplain diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..0b3f796 --- /dev/null +++ b/.gitignore @@ -0,0 +1,375 @@ +## Ignore Visual Studio temporary files, build results, and +## files generated by popular Visual Studio add-ons. +## +## Get latest from https://github.com/github/gitignore/blob/master/VisualStudio.gitignore + +# User-specific files +*.rsuser +*.suo +*.user +*.userosscache +*.sln.docstates + +# User-specific files (MonoDevelop/Xamarin Studio) +*.userprefs + +# Mono auto generated files +mono_crash.* + +# Build results +[Dd]ebug/ +[Dd]ebugPublic/ +[Rr]elease/ +[Rr]eleases/ +x64/ +x86/ +[Ww][Ii][Nn]32/ +[Aa][Rr][Mm]/ +[Aa][Rr][Mm]64/ +bld/ +[Bb]in/ +[Oo]bj/ +[Oo]ut/ +[Ll]og/ +[Ll]ogs/ + +# Visual Studio 2015/2017 cache/options directory +.vs/ +# Uncomment if you have tasks that create the project's static files in wwwroot +#wwwroot/ + +# Visual Studio 2017 auto generated files +Generated\ Files/ + +# MSTest test Results +[Tt]est[Rr]esult*/ +[Bb]uild[Ll]og.* + +# NUnit +*.VisualState.xml +TestResult.xml +nunit-*.xml + +# Build Results of an ATL Project +[Dd]ebugPS/ +[Rr]eleasePS/ +dlldata.c + +# Benchmark Results +BenchmarkDotNet.Artifacts/ + +# .NET Core +project.lock.json +project.fragment.lock.json +artifacts/ + +# ASP.NET Scaffolding +ScaffoldingReadMe.txt + +# StyleCop +StyleCopReport.xml + +# Files built by Visual Studio +*_i.c +*_p.c +*_h.h +*.ilk +*.meta +*.obj +*.iobj +*.pch +*.pdb +*.ipdb +*.pgc +*.pgd +*.rsp +*.sbr +*.tlb +*.tli +*.tlh +*.tmp +*.tmp_proj +*_wpftmp.csproj +*.log +*.vspscc +*.vssscc +.builds +*.pidb +*.svclog +*.scc +.venv + +# Chutzpah Test files +_Chutzpah* + +# Visual C++ cache files +ipch/ +*.aps +*.ncb +*.opendb +*.opensdf +*.sdf +*.cachefile +*.VC.db +*.VC.VC.opendb + +# Visual Studio profiler +*.psess +*.vsp +*.vspx +*.sap + +# Visual Studio Trace Files +*.e2e + +# TFS 2012 Local Workspace +$tf/ + +# Guidance Automation Toolkit +*.gpState + +# ReSharper is a .NET coding add-in +_ReSharper*/ +*.[Rr]e[Ss]harper +*.DotSettings.user + +# TeamCity is a build add-in +_TeamCity* + +# DotCover is a Code Coverage Tool +*.dotCover + +# AxoCover is a Code Coverage Tool +.axoCover/* +!.axoCover/settings.json + +# Coverlet is a free, cross platform Code Coverage Tool +coverage*.json +coverage*.xml +coverage*.info + +# Visual Studio code coverage results +*.coverage +*.coveragexml + +# NCrunch +_NCrunch_* +.*crunch*.local.xml +nCrunchTemp_* + +# MightyMoose +*.mm.* +AutoTest.Net/ + +# Web workbench (sass) +.sass-cache/ + +# Installshield output folder +[Ee]xpress/ + +# DocProject is a documentation generator add-in +DocProject/buildhelp/ +DocProject/Help/*.HxT +DocProject/Help/*.HxC +DocProject/Help/*.hhc +DocProject/Help/*.hhk +DocProject/Help/*.hhp +DocProject/Help/Html2 +DocProject/Help/html + +# Click-Once directory +publish/ + +# Publish Web Output +*.[Pp]ublish.xml +*.azurePubxml +# Note: Comment the next line if you want to checkin your web deploy settings, +# but database connection strings (with potential passwords) will be unencrypted +*.pubxml +*.publishproj + +# Microsoft Azure Web App publish settings. Comment the next line if you want to +# checkin your Azure Web App publish settings, but sensitive information contained +# in these scripts will be unencrypted +PublishScripts/ + +# NuGet Packages +*.nupkg +# NuGet Symbol Packages +*.snupkg +# The packages folder can be ignored because of Package Restore +**/[Pp]ackages/* +# except build/, which is used as an MSBuild target. +!**/[Pp]ackages/build/ +# Uncomment if necessary however generally it will be regenerated when needed +#!**/[Pp]ackages/repositories.config +# NuGet v3's project.json files produces more ignorable files +*.nuget.props +*.nuget.targets + +# Microsoft Azure Build Output +csx/ +*.build.csdef + +# Microsoft Azure Emulator +ecf/ +rcf/ + +# Windows Store app package directories and files +AppPackages/ +BundleArtifacts/ +Package.StoreAssociation.xml +_pkginfo.txt +*.appx +*.appxbundle +*.appxupload + +# Visual Studio cache files +# files ending in .cache can be ignored +*.[Cc]ache +# but keep track of directories ending in .cache +!?*.[Cc]ache/ + +# Others +ClientBin/ +~$* +*~ +*.dbmdl +*.dbproj.schemaview +*.jfm +*.pfx +*.publishsettings +orleans.codegen.cs + +# Including strong name files can present a security risk +# (https://github.com/github/gitignore/pull/2483#issue-259490424) +#*.snk + +# Since there are multiple workflows, uncomment next line to ignore bower_components +# (https://github.com/github/gitignore/pull/1529#issuecomment-104372622) +#bower_components/ + +# RIA/Silverlight projects +Generated_Code/ + +# Backup & report files from converting an old project file +# to a newer Visual Studio version. Backup files are not needed, +# because we have git ;-) +_UpgradeReport_Files/ +Backup*/ +UpgradeLog*.XML +UpgradeLog*.htm +ServiceFabricBackup/ +*.rptproj.bak + +# SQL Server files +*.mdf +*.ldf +*.ndf + +# Business Intelligence projects +*.rdl.data +*.bim.layout +*.bim_*.settings +*.rptproj.rsuser +*- [Bb]ackup.rdl +*- [Bb]ackup ([0-9]).rdl +*- [Bb]ackup ([0-9][0-9]).rdl + +# Microsoft Fakes +FakesAssemblies/ + +# GhostDoc plugin setting file +*.GhostDoc.xml + +# Node.js Tools for Visual Studio +.ntvs_analysis.dat +node_modules/ + +# Visual Studio 6 build log +*.plg + +# Visual Studio 6 workspace options file +*.opt + +# Visual Studio 6 auto-generated workspace file (contains which files were open etc.) +*.vbw + +# Visual Studio LightSwitch build output +**/*.HTMLClient/GeneratedArtifacts +**/*.DesktopClient/GeneratedArtifacts +**/*.DesktopClient/ModelManifest.xml +**/*.Server/GeneratedArtifacts +**/*.Server/ModelManifest.xml +_Pvt_Extensions + +# Paket dependency manager +.paket/paket.exe +paket-files/ + +# FAKE - F# Make +.fake/ + +# CodeRush personal settings +.cr/personal + +# Python Tools for Visual Studio (PTVS) +__pycache__/ +*.pyc + +# Cake - Uncomment if you are using it +# tools/** +# !tools/packages.config + +# Tabs Studio +*.tss + +# Telerik's JustMock configuration file +*.jmconfig + +# BizTalk build output +*.btp.cs +*.btm.cs +*.odx.cs +*.xsd.cs + +# OpenCover UI analysis results +OpenCover/ + +# Azure Stream Analytics local run output +ASALocalRun/ + +# MSBuild Binary and Structured Log +*.binlog + +# NVidia Nsight GPU debugger configuration file +*.nvuser + +# MFractors (Xamarin productivity tool) working folder +.mfractor/ + +# Local History for Visual Studio +.localhistory/ + +# BeatPulse healthcheck temp database +healthchecksdb + +# Backup folder for Package Reference Convert tool in Visual Studio 2017 +MigrationBackup/ + +# Ionide (cross platform F# VS Code tools) working folder +.ionide/ + +# Fody - auto-generated XML schema +FodyWeavers.xsd + +data/ +*.csv + +# Ignoruj wyuczone modele (wagi) +models/ +*.pth + +# Ignoruj wygenerowane animacje +animations/ +*.html \ No newline at end of file diff --git a/CLAUDE.md b/CLAUDE.md new file mode 100644 index 0000000..c39ad35 --- /dev/null +++ b/CLAUDE.md @@ -0,0 +1,92 @@ +# CLAUDE.md + +This file provides guidance to Claude Code (claude.ai/code) when working with code in this repository. + +## Project Overview + +**Heat Equation PINN** — A Physics-Informed Neural Network that solves the 1D time-varying heat equation with physical boundary conditions: + +``` +∂T/∂t = α ∂²T/∂x² x ∈ [0, L], t ∈ [0, T_END] +``` + +- **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=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. + +All physical and numerical parameters live in `config.py`. + +## Running + +Always activate the virtual environment first: + +```bash +source .venv/bin/activate +``` + +**PINN:** +```bash +python app.py # Train / Evaluate (L2 vs FDM) / Visualize +``` + +**FDM reference solver:** +```bash +python fdm/app.py # Solve / Heatmap / Animation / Time-series +``` + +Saved artifacts (git-ignored): `models/best_heat_pinn_model.pth`, HTML plots in `animations/` and `animations/fdm/`. + +To retrain from scratch: `rm models/best_heat_pinn_model.pth` before running option 1. + +## Dependencies + +`requirements.txt` exists. Key packages: `torch`, `numpy`, `plotly`. No `pandas` or `scikit-learn` needed. + +GPU is auto-detected (`cuda` → `mps` → `cpu`) in `engine.py:_get_device()`. + +## Architecture + +``` +config.py ← all physical + numerical parameters (edit here to change the problem) +app.py ← PINN CLI menu +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 + visualizer.py ← same 3 plot types for FDM-only output + app.py ← FDM CLI menu +``` + +### Neural Network (`model.py`) + +`HeatPINN`: 5-layer fully connected, input `(x, t)` → output `T`. + +**Output scaling** — the network predicts a dimensionless perturbation; the `forward()` applies: +``` +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. + +### 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. + +`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). + +`evaluate_model()` runs the FDM solver and downsamples its `(NX, NT)` output to the PINN prediction grid `(100, 100)` for L2 comparison. + +### 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]` +- CFL check at startup (warns, does not crash) + +### Loss Scaling Notes + +If you change `Q_VAL`, `K`, `H_CONV`, or `L` in `config.py`, the normalization in `heat_pinn_loss()` adjusts automatically. If losses diverge, check that `T_char = Q_VAL·L/K` is not near zero. diff --git a/LICENSE.txt b/LICENSE.txt new file mode 100644 index 0000000..8aa2645 --- /dev/null +++ b/LICENSE.txt @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) [year] [fullname] + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/README.md b/README.md new file mode 100644 index 0000000..9fdac78 --- /dev/null +++ b/README.md @@ -0,0 +1,2 @@ +# Projekt_NEO_PINN +# For the project to work as intended you need to download https://drive.google.com/file/d/1qZbRrY0va7xsfj1w-kk6C03Vo4LoD4ef/view?usp=drive_link and put the NEO_Curated.csv in data/ folder diff --git a/app.py b/app.py new file mode 100644 index 0000000..44bbada --- /dev/null +++ b/app.py @@ -0,0 +1,62 @@ +import sys +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("=" * 55) + + +def _ask_float(prompt, default): + val = input(prompt).strip() + try: + return float(val) + except ValueError: + return default + + +def _ask_int(prompt, default): + val = input(prompt).strip() + return int(val) if val.isdigit() else default + + +def main_menu(): + print("\nInitializing collocation points...") + data = engine.prepare_data() + print(f"Ready — device: {data['device']}\n") + + while True: + print("\n" + "-" * 30) + print(" MAIN MENU") + print("-" * 30) + print("1. Train New Model") + print("2. Evaluate Model (L2 vs analytical)") + print("3. Visualize Temperature Field") + print("0. Exit") + print("-" * 30) + + choice = input("Select an option (0-3): ").strip() + + if choice == '1': + epochs = _ask_int("Epochs (default 5000): ", 5000) + engine.train_model(data, epochs=epochs) + + elif choice == '2': + engine.evaluate_model(data) + + elif choice == '3': + engine.generate_visualization(data) + + elif choice == '0': + print("Exiting.") + sys.exit(0) + + else: + print("Invalid choice.") + + +if __name__ == "__main__": + print_header() + main_menu() diff --git a/config.py b/config.py new file mode 100644 index 0000000..32cdc6d --- /dev/null +++ b/config.py @@ -0,0 +1,20 @@ +# Fisica +ALPHA = 0.01 # diffusività termica [m²/s] +K = 1.0 # conducibilità termica [W/m·K] +L = 1.0 # lunghezza barra [m] +T0 = 20.0 # temperatura iniziale uniforme [°C] + +# Sorgente a x=0 (Neumann step) +Q_VAL = 150.0 # flusso di calore applicato [W/m²] +T_STEP = 0.2 # istante di attivazione flusso [s] + +# Convezione a x=L (Robin) +H_CONV = 10.0 # coefficiente convettivo [W/m²·K] +T_AMB = 20.0 # temperatura ambiente [°C] + +# Dominio temporale +T_END = 10.0 # fine simulazione [s] + +# Griglia FDM +NX = 100 # nodi spaziali +NT = 5000 # passi temporali (verifica CFL automatica) diff --git a/engine.py b/engine.py new file mode 100644 index 0000000..1afb783 --- /dev/null +++ b/engine.py @@ -0,0 +1,186 @@ +import os +import random +import numpy as np +import torch +import torch.optim as optim +from torch.optim.lr_scheduler import ReduceLROnPlateau + +import config +from model import HeatPINN, heat_pinn_loss +from visualizer import visualize_heat_field + +BASE_DIR = os.path.dirname(os.path.abspath(__file__)) +MODELS_DIR = os.path.join(BASE_DIR, 'models') +MODEL_SAVE_PATH = os.path.join(MODELS_DIR, 'best_heat_pinn_model.pth') + + +def set_seed(seed=42): + random.seed(seed) + np.random.seed(seed) + torch.manual_seed(seed) + if torch.cuda.is_available(): + torch.cuda.manual_seed_all(seed) + + +def _get_device(): + if torch.cuda.is_available(): + try: + t = torch.zeros(2, 2).cuda() + torch.mm(t, t) + return torch.device('cuda') + except RuntimeError: + pass + if torch.backends.mps.is_available(): + return torch.device('mps') + return torch.device('cpu') + + +def prepare_data(N_f=4000, N_ic=400, N_bc=400): + set_seed(42) + device = _get_device() + + # Uniform collocation points + 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) + 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_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_ic = torch.rand(N_ic, device=device) * config.L + t_bc = torch.rand(N_bc, device=device) * config.T_END + + return {'device': device, 'x_f': x_f, 't_f': t_f, 'x_ic': x_ic, 't_bc': t_bc} + + +def train_model(data, epochs=5000, patience=100): + device = data['device'] + model = HeatPINN().to(device) + optimizer = optim.Adam(model.parameters(), lr=1e-3) + scheduler = ReduceLROnPlateau(optimizer, mode='min', factor=0.5, patience=30, min_lr=1e-6) + + os.makedirs(MODELS_DIR, exist_ok=True) + best_loss = float('inf') + patience_counter = 0 + + print(f"\n--- Heat PINN Training (Adam) on {device} ---") + for epoch in range(epochs): + model.train() + optimizer.zero_grad() + loss, L_pde, L_ic, L_bc = heat_pinn_loss( + model, data['x_f'], data['t_f'], data['x_ic'], data['t_bc'] + ) + loss.backward() + optimizer.step() + scheduler.step(loss.item()) + + if loss.item() < best_loss - 1e-7: + best_loss = loss.item() + patience_counter = 0 + torch.save({'state_dict': model.state_dict()}, MODEL_SAVE_PATH) + else: + patience_counter += 1 + + if patience_counter >= patience: + print(f"Early stopping at epoch {epoch + 1}") + break + + if (epoch + 1) % 100 == 0: + print( + f"Epoch [{epoch+1}/{epochs}] | Loss: {loss.item():.6f} " + f"| PDE: {L_pde.item():.6f} | IC: {L_ic.item():.6f} | BC: {L_bc.item():.6f}" + ) + + # L-BFGS fine-tuning phase (standard PINN practice for convergence to better minima) + print("\n--- L-BFGS fine-tuning ---") + ckpt = torch.load(MODEL_SAVE_PATH, map_location=device) + model.load_state_dict(ckpt['state_dict']) + + lbfgs = optim.LBFGS(model.parameters(), lr=0.1, max_iter=50, + history_size=50, tolerance_grad=1e-7, line_search_fn='strong_wolfe') + + _last = {} + + for step in range(20): + def closure(): + lbfgs.zero_grad() + loss, L_pde, L_ic, L_bc = heat_pinn_loss( + model, data['x_f'], data['t_f'], data['x_ic'], data['t_bc'] + ) + loss.backward() + _last['loss'] = loss.item() + _last['pde'] = L_pde.item() + _last['ic'] = L_ic.item() + _last['bc'] = L_bc.item() + return loss + + lbfgs.step(closure) + if _last['loss'] < best_loss: + best_loss = _last['loss'] + torch.save({'state_dict': model.state_dict()}, MODEL_SAVE_PATH) + if (step + 1) % 5 == 0: + print(f"L-BFGS step {step+1}/20 | Loss: {_last['loss']:.6f} " + f"| PDE: {_last['pde']:.6f} | IC: {_last['ic']:.6f} | BC: {_last['bc']:.6f}") + + print("Training complete! Model saved.") + + +def _load_model(device): + if not os.path.exists(MODEL_SAVE_PATH): + print("Error: model not found. Train the model first.") + return None + ckpt = torch.load(MODEL_SAVE_PATH, map_location=device) + model = HeatPINN().to(device) + model.load_state_dict(ckpt['state_dict']) + model.eval() + return model + + +def _predict_grid(model, device, nx=100, nt=100): + x_vals = torch.linspace(0, config.L, nx, device=device) + t_vals = torch.linspace(0, config.T_END, nt, device=device) + xx, tt = torch.meshgrid(x_vals, t_vals, indexing='ij') + inp = torch.stack([xx.reshape(-1), tt.reshape(-1)], dim=1) + with torch.no_grad(): + T_pred = model(inp).reshape(nx, nt).cpu().numpy() + return T_pred, x_vals.cpu().numpy(), t_vals.cpu().numpy() + + +def evaluate_model(data): + model = _load_model(data['device']) + if model is None: + return + T_pred, x_vals, t_vals = _predict_grid(model, data['device']) + + # FDM reference + 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) + nx, nt = T_pred.shape + t_indices = np.linspace(0, T_fdm.shape[1] - 1, nt, dtype=int) + T_fdm_ds = T_fdm[:, 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)) + + print(f"\n--- Evaluation vs FDM ---") + print(f"Relative L2 error vs FDM: {l2_rel:.6f}") + print(f"Max absolute error: {max_err:.6f} °C\n") + + +def generate_visualization(data): + model = _load_model(data['device']) + if model is None: + return + T_pred, x_vals, t_vals = _predict_grid(model, data['device']) + from fdm.solver import solve as fdm_solve + T_fdm, _, _ = fdm_solve() + visualize_heat_field(T_pred, x_vals, t_vals, T_fdm) diff --git a/fdm/__init__.py b/fdm/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/fdm/app.py b/fdm/app.py new file mode 100644 index 0000000..b549db5 --- /dev/null +++ b/fdm/app.py @@ -0,0 +1,58 @@ +import sys +import os + +sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) + +from fdm.solver import solve +from fdm.visualizer import visualize_fdm + + +def print_header(): + print("=" * 38) + print(" Heat Equation — FDM Solver") + print(" ∂T/∂t = α ∂²T/∂x²") + print("=" * 38) + + +def main_menu(): + print("\nInitializing solver...") + print("Ready.\n") + + T, x_vals, t_vals = None, None, None + + while True: + 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("0. Esci") + print("-" * 30) + + choice = input("Select an option (0-4): ").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"): + 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.") + sys.exit(0) + + else: + print("Scelta non valida.") + + +if __name__ == "__main__": + print_header() + main_menu() diff --git a/fdm/solver.py b/fdm/solver.py new file mode 100644 index 0000000..f703e17 --- /dev/null +++ b/fdm/solver.py @@ -0,0 +1,103 @@ +""" +FTCS (Forward-Time Centered-Space) explicit finite difference solver +for the 1D heat equation: + + dT/dt = alpha * d²T/dx² + +Boundary conditions: + - x=0 (Neumann): heat flux step Q(t) applied via ghost cell + - x=L (Robin): convective boundary condition + +Returns T_matrix of shape (NX, NT). +""" + +import sys +import os +import numpy as np + +# Allow importing config from the project root +sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) +import config + + +def solve(): + """Run the FTCS solver and return (T_matrix, x_vals, t_vals). + + Returns + ------- + T_matrix : np.ndarray, shape (NX, NT) + Temperature at each spatial node (row) and time step (column). + T_matrix[i, n] = T(x_i, t_n). + x_vals : np.ndarray, shape (NX,) + Spatial node positions from 0 to L. + t_vals : np.ndarray, shape (NT,) + Time values from 0 to T_END. + """ + # ----------------------------------------------------------------- + # Parameters from config + # ----------------------------------------------------------------- + alpha = config.ALPHA + k = config.K + L = config.L + T0 = config.T0 + Q_val = config.Q_VAL + t_step = config.T_STEP + h = config.H_CONV + T_amb = config.T_AMB + T_end = config.T_END + NX = config.NX + NT = config.NT + + # ----------------------------------------------------------------- + # Grid + # ----------------------------------------------------------------- + x_vals = np.linspace(0.0, L, NX) + t_vals = np.linspace(0.0, T_end, NT) + dx = x_vals[1] - x_vals[0] + dt = t_vals[1] - t_vals[0] + + # ----------------------------------------------------------------- + # Stability check (CFL for explicit heat equation: r <= 0.5) + # ----------------------------------------------------------------- + r = alpha * dt / dx**2 + if dt > dx**2 / (2.0 * alpha): + print( + f"[FDM WARNING] Stability condition violated: " + f"dt={dt:.6g} > dx²/(2*alpha)={dx**2/(2.0*alpha):.6g} (r={r:.4f} > 0.5). " + "Solution may diverge." + ) + + # ----------------------------------------------------------------- + # Allocate output matrix and set initial condition + # ----------------------------------------------------------------- + T_matrix = np.zeros((NX, NT), dtype=np.float64) + T_matrix[:, 0] = T0 # uniform IC + + # Working array for the current time level + T_cur = np.full(NX, T0, dtype=np.float64) + + # ----------------------------------------------------------------- + # Time integration + # ----------------------------------------------------------------- + for n in range(NT - 1): + t_now = t_vals[n] + + # --- Neumann BC at x=0: ghost cell --- + # Q(t) = Q_val if t >= t_step else 0 + Q = Q_val if t_now >= t_step else 0.0 + T_ghost = T_cur[1] + 2.0 * dx * Q / k # ghost node T[-1] + + # --- Interior FTCS update (indices 1 .. NX-2) --- + T_new = T_cur.copy() + T_new[1:-1] = T_cur[1:-1] + r * (T_cur[2:] - 2.0 * T_cur[1:-1] + T_cur[:-2]) + + # --- Apply Neumann BC at i=0 using ghost cell --- + T_new[0] = T_cur[0] + r * (T_cur[1] - 2.0 * T_cur[0] + T_ghost) + + # --- Apply Robin BC at x=L (explicit) --- + T_new[-1] = (T_cur[-2] + dx * h / k * T_amb) / (1.0 + dx * h / k) + + T_cur = T_new + T_matrix[:, n + 1] = T_cur + + return T_matrix, x_vals, t_vals diff --git a/fdm/visualizer.py b/fdm/visualizer.py new file mode 100644 index 0000000..37e2334 --- /dev/null +++ b/fdm/visualizer.py @@ -0,0 +1,227 @@ +import sys +import os +import numpy as np +import plotly.graph_objects as go + +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): + """Produce three HTML visualisations for the FDM solver output. + + Parameters + ---------- + T_matrix : np.ndarray, shape (NX, NT) + Temperature field: T_matrix[i, n] = T(x_i, t_n). + x_vals : np.ndarray, shape (NX,) + Spatial node positions. + t_vals : np.ndarray, shape (NT,) + Time values. + """ + os.makedirs(FDM_ANIM_DIR, exist_ok=True) + + # ------------------------------------------------------------------ + # 1. Heatmap + # ------------------------------------------------------------------ + 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 + + 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)) + + fig_heatmap = go.Figure(go.Heatmap( + z=z_data, + x=x_vals, + y=t_vals, + colorscale='RdBu_r', + zmin=zmin_sym, + zmax=zmax_sym, + colorbar=dict(title='T [°C]'), + )) + fig_heatmap.update_layout( + title='FDM — Heat Equation T(x,t)', + xaxis_title='x [m]', + yaxis_title='t [s]', + height=500, + ) + + heatmap_path = _next_path(FDM_ANIM_DIR, 'heatmap', '.html') + fig_heatmap.write_html(heatmap_path) + 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 + + vline_shapes = [ + dict(type='line', x0=0, x1=0, y0=y_min, y1=y_max, + line=dict(color='grey', dash='dash', width=1)), + dict(type='line', x0=L, x1=L, y0=y_min, y1=y_max, + line=dict(color='grey', dash='dash', width=1)), + ] + + frames = [] + for idx in frame_indices: + frames.append(go.Frame( + data=[go.Scatter( + x=x_vals, + y=T_matrix[:, idx], + mode='lines', + line=dict(color='royalblue', width=2), + name='FDM', + )], + name=str(idx), + layout=go.Layout( + title_text=f'FDM | t = {t_vals[idx]:.3f} s', + ), + )) + + fig_anim = go.Figure( + data=[go.Scatter( + x=x_vals, + y=T_matrix[:, frame_indices[0]], + mode='lines', + line=dict(color='royalblue', width=2), + name='FDM', + )], + layout=go.Layout( + title=f'FDM | t = {t_vals[frame_indices[0]]:.3f} s', + xaxis=dict(title='x [m]', range=[-0.02 * L, 1.02 * L]), + yaxis=dict(title='T [°C]', range=[y_min, y_max]), + shapes=vline_shapes, + updatemenus=[dict( + type='buttons', + showactive=False, + y=1.15, + x=0.5, + xanchor='center', + buttons=[ + dict( + label='▶ Play', + method='animate', + args=[None, dict( + frame=dict(duration=40, redraw=False), + 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=False), + )], + 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)), + )], + ), + frames=frames, + ) + + anim_path = _next_path(FDM_ANIM_DIR, 'animation', '.html') + fig_anim.write_html(anim_path) + print(f"Animation saved → {anim_path}") + + # ------------------------------------------------------------------ + # 3. Time evolution at fixed spatial points + # ------------------------------------------------------------------ + fixed_fractions = [0.0, 0.25, 0.5, 0.75, 1.0] + colors = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd'] + + fig_ts = go.Figure() + + for frac, color in zip(fixed_fractions, colors): + x_target = frac * L + idx_x = int(np.argmin(np.abs(x_vals - x_target))) + label = f'x = {x_vals[idx_x]:.2f} m' + fig_ts.add_trace(go.Scatter( + x=t_vals, + y=T_matrix[idx_x, :], + mode='lines', + name=label, + line=dict(color=color, width=2), + )) + + # "Heat ON" vertical dashed line at t = T_STEP + t_step = config.T_STEP + t_min_val = float(np.min(t_vals)) + t_max_val = float(np.max(t_vals)) + T_min_val = float(np.min(T_matrix)) - 1.0 + T_max_val = float(np.max(T_matrix)) + 1.0 + + fig_ts.add_shape( + type='line', + x0=t_step, x1=t_step, + y0=T_min_val, y1=T_max_val, + line=dict(color='red', dash='dash', width=1.5), + ) + fig_ts.add_annotation( + x=t_step, + y=T_max_val, + text='Heat ON', + showarrow=False, + yanchor='top', + font=dict(color='red', size=11), + ) + + fig_ts.update_layout( + title='FDM — Temperature evolution at fixed points', + xaxis=dict(title='t [s]', range=[t_min_val, t_max_val]), + yaxis=dict(title='T [°C]', range=[T_min_val, T_max_val]), + legend=dict(x=0.01, y=0.99), + height=480, + ) + + ts_path = _next_path(FDM_ANIM_DIR, 'time_series', '.html') + fig_ts.write_html(ts_path) + print(f"Time-series saved → {ts_path}") diff --git a/model.py b/model.py new file mode 100644 index 0000000..6600c19 --- /dev/null +++ b/model.py @@ -0,0 +1,61 @@ +import torch +import torch.nn as nn +import config + + +class HeatPINN(nn.Module): + def __init__(self): + super().__init__() + self.net = nn.Sequential( + nn.Linear(2, 128), nn.Tanh(), + nn.Linear(128, 128), nn.Tanh(), + nn.Linear(128, 128), nn.Tanh(), + nn.Linear(128, 128), nn.Tanh(), + nn.Linear(128, 1), + ) + + def forward(self, x): + # Output scaled to physical range: T_AMB + (Q*L/K) * net + # net learns dimensionless perturbation in [0,1] range + T_scale = config.T_AMB + (config.Q_VAL * config.L / config.K) * self.net(x) + return T_scale + + +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): + # Characteristic scales for normalization + 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) + 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] + pde_scale = (T_char / config.T_END) ** 2 + 1e-8 + L_pde = ((dT_dt - config.ALPHA * d2T_dx2) ** 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 + 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 + + # 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) + T_right = model(torch.stack([x_right, t_bc.detach()], dim=1)) + dT_dx_right = torch.autograd.grad(T_right.sum(), x_right, create_graph=True)[0] + L_bc_right = ((dT_dx_right + (config.H_CONV / config.K) * (T_right.squeeze() - config.T_AMB)) ** 2).mean() / grad_char + + L_bc = L_bc_left + L_bc_right + total = w_pde * L_pde + w_ic * L_ic + w_bc * L_bc + return total, L_pde, L_ic, L_bc diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..687263c --- /dev/null +++ b/requirements.txt @@ -0,0 +1,5 @@ +torch>=2.0.0 +pandas>=2.0.0 +numpy>=1.24.0 +scikit-learn>=1.3.0 +plotly>=5.15.0 diff --git a/visualizer.py b/visualizer.py new file mode 100644 index 0000000..be21687 --- /dev/null +++ b/visualizer.py @@ -0,0 +1,153 @@ +import os +import numpy as np +import plotly.graph_objects as go +from plotly.subplots import make_subplots +import config + +BASE_DIR = os.path.dirname(os.path.abspath(__file__)) +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)) + nt_pred = len(t_vals) + 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) + + # --- Static heatmap: PINN vs FDM --- + fig_map = make_subplots( + rows=1, cols=2, + subplot_titles=["PINN Prediction T(x,t)", "FDM Reference T(x,t)"], + shared_yaxes=True, + ) + colorscale = 'Hot_r' + zmin = float(min(np.min(T_pred), np.min(T_fdm_ds))) + zmax = float(max(np.max(T_pred), np.max(T_fdm_ds))) + + fig_map.add_trace(go.Heatmap( + z=T_pred.T, x=x_vals, y=t_vals, + colorscale=colorscale, zmin=zmin, zmax=zmax, + colorbar=dict(x=0.46, title='T [°C]'), + ), row=1, col=1) + + fig_map.add_trace(go.Heatmap( + z=T_fdm_ds.T, x=x_vals, y=t_vals, + colorscale=colorscale, zmin=zmin, zmax=zmax, + colorbar=dict(x=1.01, title='T [°C]'), + ), row=1, col=2) + + fig_map.update_xaxes(title_text='x') + fig_map.update_yaxes(title_text='t', row=1, col=1) + fig_map.update_layout(title_text='Heat Equation PINN vs FDM', height=450) + + map_path = _next_path('heatmap', '.html') + fig_map.write_html(map_path) + print(f"Heatmap saved → {map_path}") + + # --- Animated profile T(x) evolving in time --- + n_frames = len(t_vals) + frames = [] + for i in range(n_frames): + frames.append(go.Frame( + data=[ + go.Scatter(x=x_vals, y=T_pred[:, i], mode='lines', + line=dict(color='royalblue', width=2), name='PINN'), + go.Scatter(x=x_vals, y=T_fdm_ds[:, i], mode='lines', + line=dict(color='firebrick', width=2, dash='dash'), name='FDM'), + ], + name=str(i), + layout=go.Layout(title_text=f'Heat Equation PINN vs FDM | t = {t_vals[i]:.3f}'), + )) + + fig_anim = go.Figure( + data=frames[0].data, + layout=go.Layout( + title=f'Heat Equation PINN vs FDM | t = {t_vals[0]:.3f}', + xaxis=dict(title='x [m]', range=[-0.02, config.L + 0.02]), + yaxis=dict(title='T [°C]', range=[zmin - 1, zmax + 1]), + legend=dict(x=0.75, y=0.95), + updatemenus=[dict( + type='buttons', showactive=False, y=1.15, x=0.5, xanchor='center', + buttons=[ + dict(label='▶ Play', method='animate', + args=[None, dict(frame=dict(duration=40, redraw=False), + 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(i)], + dict(mode='immediate', frame=dict(duration=0, redraw=False))], + label=f'{t_vals[i]:.2f}') for i in range(n_frames)], + transition=dict(duration=0), x=0.05, y=0, len=0.9, + currentvalue=dict(prefix='t = ', font=dict(size=14)), + )], + ), + frames=frames, + ) + + anim_path = _next_path('heat_animation', '.html') + fig_anim.write_html(anim_path) + print(f"Animation saved → {anim_path}") + + # --- Time-series comparison at fixed spatial points --- + # Spatial indices for x=0, x=L/2, x=L + nx = len(x_vals) + idx_x0 = 0 + idx_xmid = nx // 2 + idx_xL = nx - 1 + + points = [ + (idx_x0, 'x=0', 'blue'), + (idx_xmid, 'x=L/2', 'green'), + (idx_xL, 'x=L', 'red'), + ] + + fig_ts = go.Figure() + + for idx, label, color in points: + fig_ts.add_trace(go.Scatter( + x=t_vals, y=T_pred[idx, :], + mode='lines', + line=dict(color=color, width=2, dash='solid'), + name=f'PINN {label}', + )) + fig_ts.add_trace(go.Scatter( + x=t_vals, y=T_fdm_ds[idx, :], + mode='lines', + line=dict(color=color, width=2, dash='dash'), + name=f'FDM {label}', + )) + + # Vertical dashed line at T_STEP ("Heat ON") + fig_ts.add_vline( + x=config.T_STEP, + line=dict(color='red', dash='dash', width=1.5), + annotation_text='Heat ON', + annotation_position='top right', + ) + + fig_ts.update_layout( + title='Heat Equation PINN vs FDM — Time Series at Fixed Points', + xaxis_title='t', + yaxis_title='T(x,t)', + legend=dict(x=0.01, y=0.99), + height=500, + ) + + comparison_path = _next_path('comparison', '.html') + fig_ts.write_html(comparison_path) + print(f"Time-series saved → {comparison_path}") + + +def _next_path(prefix, ext): + i = 1 + while True: + path = os.path.join(ANIMATIONS_DIR, f'{prefix}_{i:03d}{ext}') + if not os.path.exists(path): + return path + i += 1