From 0f4438cc7fbf7c03645765d221ebda34339a232f Mon Sep 17 00:00:00 2001 From: Davide Grilli Date: Wed, 13 May 2026 21:21:53 +0200 Subject: [PATCH] Commit iniziale --- .gitattributes | 63 ++++++++ .gitignore | 375 ++++++++++++++++++++++++++++++++++++++++++++++ CLAUDE.md | 92 ++++++++++++ LICENSE.txt | 21 +++ README.md | 2 + app.py | 62 ++++++++ config.py | 20 +++ engine.py | 186 +++++++++++++++++++++++ fdm/__init__.py | 0 fdm/app.py | 58 +++++++ fdm/solver.py | 103 +++++++++++++ fdm/visualizer.py | 227 ++++++++++++++++++++++++++++ model.py | 61 ++++++++ requirements.txt | 5 + visualizer.py | 153 +++++++++++++++++++ 15 files changed, 1428 insertions(+) create mode 100644 .gitattributes create mode 100644 .gitignore create mode 100644 CLAUDE.md create mode 100644 LICENSE.txt create mode 100644 README.md create mode 100644 app.py create mode 100644 config.py create mode 100644 engine.py create mode 100644 fdm/__init__.py create mode 100644 fdm/app.py create mode 100644 fdm/solver.py create mode 100644 fdm/visualizer.py create mode 100644 model.py create mode 100644 requirements.txt create mode 100644 visualizer.py 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