diff --git a/DOCUMENTATION_STYLE.md b/DOCUMENTATION_STYLE.md index e1367bd3e..cfe9fc175 100644 --- a/DOCUMENTATION_STYLE.md +++ b/DOCUMENTATION_STYLE.md @@ -340,8 +340,6 @@ Prefer a more specific term whenever one exists. | AT1 bound treatment | projected conjugate gradients, then projected CG | | Explicit route | explicit dynamics with the integrator named | | Quasi-static route | implicit equilibrium solve within the staggered scheme | -| Learned proposal | learned damage proposal followed by the documented checks | -| Learned replacement | audited learned damage replacement with classical fallback | | Subcycling | damage subcycling with interval `$N$` | | Fracture toughness | `$G_c$`, with units defined | | Regularisation length | `$\ell_0$`, with units defined | @@ -384,7 +382,7 @@ When comparing performance, report: Do not describe one backend as universally fastest. -## 10. Differentiability and learned models +## 10. Differentiability Avoid unqualified statements that PhAST differentiates every fracture calculation. Prefer: @@ -393,18 +391,6 @@ calculation. Prefer: > updates, bounds, active-set changes, remeshing decisions, and other > non-smooth operations require case-specific interpretation. -For learned damage updates, state who remains responsible for the physical -solution: - -- A learned proposal supplies an initial damage estimate. -- PhAST checks bounds, irreversibility, phase-field boundary conditions, and - the projected residual where the selected route provides those checks. -- The classical damage solver remains available as the fallback. -- No trained model is implied unless a checkpoint is distributed and linked. - -Do not call a learned model faster without measuring its complete inference, -checking, rejection, and fallback cost against the matched classical solve. - ## 11. Commands and configurations - Copy commands from the current CLI help. diff --git a/README.md b/README.md index ff950dd86..e5065d0ca 100644 --- a/README.md +++ b/README.md @@ -223,7 +223,6 @@ Browse the full [example gallery](docs/example-gallery.md) for the complete list | Post-process simulation data | `phast.load_result(path)` | [Public API Reference](docs/user_guide/public_api_reference.md) | | Review supported physics | Capability matrix | [Capability Matrix](docs/user_guide/capability_matrix.md) | | Learn step-by-step setup | Tutorial notebook | [Problem Setup Walkthrough](docs/tutorial/problem_setup_walkthrough.ipynb) | -| Add an audited learned damage model | Predictor plug-in protocol | [Modular FEM and Learned Damage](docs/tutorial/03_modular_fem_and_learned_damage.md) | | Learn elementwise `E(x)` and `Gc(x)` fields | Script-contract teaching example | [Heterogeneous Material Fields](docs/tutorial/05_heterogeneous_material_fields.md) | | Diagnose failed runs | Troubleshooting guide | [Troubleshooting](docs/troubleshooting.md) | diff --git a/configs/REFERENCE.yaml b/configs/REFERENCE.yaml index 6594a436c..2239525bc 100644 --- a/configs/REFERENCE.yaml +++ b/configs/REFERENCE.yaml @@ -386,14 +386,6 @@ solver: # softmax_H_beta: null # float or null # H_update_method: hard_max # str, enum: hard_max|softmax|smooth_max|log_smooth|custom_subgrad # enable_damage: true # bool - # damage_update: classical # str, enum: classical|learned_proposal|learned_replacement - # damage_predictor: null # str or null, module:factory used by learned routes - # damage_checkpoint: null # str or null, interpreted by the predictor factory - # damage_predictor_options: {} # mapping, architecture-specific adapter options - # damage_residual_rtol: 0.001 # float, range: [0.0, +inf] - # damage_residual_atol: 1e-08 # float, range: [0.0, +inf] - # damage_bound_tolerance: 1e-08 # float, range: [0.0, +inf] - # damage_fallback: true # bool, use classical damage solve after rejection # fail_on_mechanics_nonconvergence: true # bool # adaptive_dt: false # bool # adaptive_dt_d_threshold: 0.01 # float diff --git a/configs/phast.schema.json b/configs/phast.schema.json index 8b853691f..6a6dac9d3 100644 --- a/configs/phast.schema.json +++ b/configs/phast.schema.json @@ -152,7 +152,14 @@ "enum": [ "standard", "cubic", - "rational" + "rational", + "rational_at2" + ] + }, + "stress_degradation": { + "enum": [ + "split", + "full" ] }, "driving_force": { @@ -308,7 +315,19 @@ "enum": [ "standard", "cubic", - "rational" + "rational", + "rational_at2" + ] + }, + "stress_degradation": { + "type": [ + "string", + "null" + ], + "default": null, + "enum": [ + "split", + "full" ] }, "kinematics": { @@ -622,53 +641,6 @@ "type": "boolean", "default": true }, - "damage_update": { - "type": "string", - "default": "classical", - "enum": [ - "classical", - "learned_proposal", - "learned_replacement" - ] - }, - "damage_predictor": { - "type": [ - "string", - "null" - ], - "default": null - }, - "damage_checkpoint": { - "type": [ - "string", - "null" - ], - "default": null - }, - "damage_predictor_options": { - "type": "object", - "additionalProperties": true, - "default": {} - }, - "damage_residual_rtol": { - "type": "number", - "default": 0.001, - "minimum": 0.0 - }, - "damage_residual_atol": { - "type": "number", - "default": 1e-08, - "minimum": 0.0 - }, - "damage_bound_tolerance": { - "type": "number", - "default": 1e-08, - "minimum": 0.0 - }, - "damage_fallback": { - "type": "boolean", - "default": true - }, "fail_on_mechanics_nonconvergence": { "type": "boolean", "default": true diff --git a/docs/index.md b/docs/index.md index 8d5081f78..7d0c12d75 100644 --- a/docs/index.md +++ b/docs/index.md @@ -148,7 +148,7 @@ tutorial/notebook_dynamic_sent tutorial/notebook_setup tutorial/notebook_mesh_resolution tutorial/notebook_retained_results -tutorial/03_modular_fem_and_learned_damage +tutorial/03_modular_fem tutorial/04_exploration_experiments tutorial/05_heterogeneous_material_fields tutorial/06_first_research_study @@ -176,7 +176,6 @@ user_guide/setup_problems user_guide/python_api user_guide/yaml_workflow user_guide/results_visualization -user_guide/learned_damage performance-reproducibility paper-reproduction ``` diff --git a/docs/tutorial/03_modular_fem.md b/docs/tutorial/03_modular_fem.md new file mode 100644 index 000000000..9277b56cd --- /dev/null +++ b/docs/tutorial/03_modular_fem.md @@ -0,0 +1,94 @@ +# Modular fracture problems + +This tutorial describes how the main components of a phase-field fracture +problem fit together in PhAST without assuming prior knowledge of the +codebase. + +| Lesson item | Scope | +|---|---| +| Prerequisites | Complete the installation checks and read the phase-field primer. | +| Route | Quasi-static staggered phase-field fracture through a checked-in YAML configuration. | +| Hardware | CPU is sufficient for configuration inspection and the bounded example route. | +| Verification | Inspect the selected route, result manifest, and convergence history. | + +## 1. The finite-element problem + +A runnable problem combines: + +1. A mesh defining geometry, nodes, elements, and physical boundary groups. +2. A material defining elasticity, fracture resistance, regularization length, + density, crack-density model, degradation law, and energy split. +3. Initial conditions defining displacement, velocity, damage, and any initial + crack. +4. Boundary conditions constraining mechanical or phase-field degrees of + freedom. +5. Loading prescribing the time-dependent displacement or force history. +6. Solver settings choosing integration, linear algebra, and the damage + subproblem. +7. Output settings selecting fields and sampling frequency. + +Keep reusable configurations under `configs/`. Schema, manifest, contract, and +template files are machine-readable specifications rather than simulations. + +## 2. Distinguish the fracture-model choices + +- `pf_model: AT1` or `AT2` selects the crack-density functional. +- `energy_split` selects the elastic-energy decomposition. Implementations + include `spectral`, `spectral_stress`, + `spectral_plane_stress_condensed`, `isotropic`, `amor`, and `star_convex`. +- `degradation_type` selects the stiffness-degradation family. + +The supported coupled AT1/AT2 route uses `degradation_type: standard`. +PF-CZM is experimental and constructs its rational law from its own +parameters. Although the material layer contains additional degradation +functions, PhAST does not present cubic or rational degradation as supported +coupled AT1/AT2 damage solves. + +## 3. Begin from a complete problem + +Mesh groups, loading controls, and output requests are geometry-specific. Use +the checked-in Miehe tension configuration rather than a fragment containing +placeholder paths: + +```bash +python -m phast explain-config examples/quasistatic/miehe_tension/config.yaml +python -m phast run examples/quasistatic/miehe_tension/config.yaml --validate-only +python -m phast run examples/quasistatic/miehe_tension/config.yaml \ + --num_steps 5 --output_dir runs/miehe_tension_five_step +``` + +The five-step calculation is a bounded workflow exercise, not a crack-growth +validation result. Read the example README before undertaking the full +calculation; it records the retained runtime and evidence boundary. + +## 4. How automatic route selection works + +Automatic device and sparse-backend settings select a deterministic route from +the backends available in the current environment. This is a portability +policy, not a universal fastest-solver guarantee. Performance depends on +problem size, sparsity, hardware, installed libraries, and transfer costs. + +At startup and completion, PhAST reports the device, precision, element type, +mechanics solver, configured and resolved sparse backend, fracture model, +energy split, history update, damage-bound treatment, and preconditioner. +Preserve this report with benchmark results. + +```bash +python -m phast doctor +``` + +## 5. A defensible student workflow + +1. Run an existing example without changing its formulation. +2. Modify one component at a time: mesh, material, loading, or energy split. +3. Examine the route report and retained displacement, damage, and history. +4. Compare each modified run with the unmodified reference calculation. + +Students and developers are invited to propose tutorials and examples. If an +installation step, configuration field, or solver message is unclear, open a +GitHub issue with the configuration, environment report, and minimal +reproduction. + +Continue with [Heterogeneous material fields](05_heterogeneous_material_fields.md) +for an element-ordered material-field teaching problem, or return to the +[example gallery](../example-gallery.md) for complete workflows. diff --git a/docs/tutorial/03_modular_fem_and_learned_damage.md b/docs/tutorial/03_modular_fem_and_learned_damage.md deleted file mode 100644 index d2e28bf2a..000000000 --- a/docs/tutorial/03_modular_fem_and_learned_damage.md +++ /dev/null @@ -1,164 +0,0 @@ -# Modular fracture problems and learned damage updates - -This tutorial describes how the main components of a phase-field fracture -problem fit together in PhAST. It also introduces the experimental -learned-damage plug-in without assuming prior knowledge of the codebase. - -| Lesson item | Scope | -|---|---| -| Prerequisites | Complete the installation checks and read the phase-field primer. | -| Classical route | Quasi-static staggered phase-field fracture through a checked-in YAML configuration. | -| Experimental route | A user-supplied learned proposal or replacement for the damage subproblem only. | -| Hardware | CPU is sufficient for configuration inspection and the bounded example route. | -| Verification | Inspect the selected route, result manifest, convergence history, and learned acceptance/fallback telemetry. | - -## 1. The finite-element problem - -A runnable problem combines: - -1. A mesh defining geometry, nodes, elements, and physical boundary groups. -2. A material defining elasticity, fracture resistance, regularization length, - density, crack-density model, degradation law, and energy split. -3. Initial conditions defining displacement, velocity, damage, and any initial - crack. -4. Boundary conditions constraining mechanical or phase-field degrees of - freedom. -5. Loading prescribing the time-dependent displacement or force history. -6. Solver settings choosing integration, linear algebra, and damage update. -7. Output settings selecting fields and sampling frequency. - -Keep reusable configurations under `configs/`. Schema, manifest, contract, and -template files are machine-readable specifications rather than simulations. - -## 2. Distinguish the fracture-model choices - -- `pf_model: AT1` or `AT2` selects the crack-density functional. -- `energy_split` selects the elastic-energy decomposition. Implementations - include `spectral`, `spectral_stress`, - `spectral_plane_stress_condensed`, `isotropic`, `amor`, and `star_convex`. -- `degradation_type` selects the stiffness-degradation family. - -The supported coupled AT1/AT2 route uses `degradation_type: standard`. -PF-CZM is experimental and constructs its rational law from its own -parameters. Although the material layer contains additional degradation -functions, PhAST does not present cubic or rational degradation as supported -coupled AT1/AT2 damage solves. - -## 3. Begin from a complete classical problem - -Mesh groups, loading controls, and output requests are geometry-specific. Use -the checked-in Miehe tension configuration rather than a fragment containing -placeholder paths: - -```bash -python -m phast explain-config examples/quasistatic/miehe_tension/config.yaml -python -m phast run examples/quasistatic/miehe_tension/config.yaml --validate-only -python -m phast run examples/quasistatic/miehe_tension/config.yaml \ - --num_steps 5 --output_dir runs/miehe_tension_five_step -``` - -The five-step calculation is a bounded workflow exercise, not a crack-growth -validation result. Read the example README before undertaking the full -calculation; it records the retained runtime and evidence boundary. - -## 4. How automatic route selection works - -Automatic device and sparse-backend settings select a deterministic route from -the backends available in the current environment. This is a portability -policy, not a universal fastest-solver guarantee. Performance depends on -problem size, sparsity, hardware, installed libraries, and transfer costs. - -At startup and completion, PhAST reports the device, configured and resolved -sparse backend, damage-update mode, predictor, and learned-update -acceptance/fallback counts. Preserve this report with benchmark results. - -```bash -python -m phast doctor -``` - -## 5. Use a learned model as a proposal - -PhAST does not distribute a generally applicable trained damage checkpoint. -The following block documents the plug-in contract and is not independently -runnable until the user supplies a compatible factory and checkpoint. - -```yaml -solver: - damage_update: learned_proposal - damage_predictor: examples.learned_damage.predictor_plugin:create_predictor - damage_checkpoint: path/to/model.pt - damage_predictor_options: - representation: damage_increment -``` - -The predictor receives a `DamageStepContext` containing coordinates, -connectivity, displacement, velocity, element and nodal history, previous -damage, material metadata, formulation choices, time, and load factor. The -reference TorchScript adapter uses: - -```text -[x, y, history field, previous damage, displacement x, displacement y] -``` - -This six-column matrix is only the reference adapter. The public protocol does -not impose one architecture or tensor layout. A plug-in can access raw -coordinates and connectivity, call `context.graph_edge_index()`, construct -node-element incidence operators, retain temporal state, project fields to a -regular grid, and apply checkpoint-specific normalization. This accommodates -graph, coordinate-conditioned, temporal, grid-operator, and reduced-basis -research without placing Paper-specific models in the solver core. - -A `DamagePrediction` declares `representation="damage"` or -`representation="damage_increment"`. The latter is added to the previous -accepted damage before projection and audit. - -PhAST projects the proposal onto damage bounds, irreversibility, and prescribed -phase-field values, then uses it as the initial guess for the classical damage -solve. The converged classical update remains authoritative. - -## 6. Experimental learned replacement - -This mode still replaces only the damage-subproblem update. Mechanics, history -construction, constraints, audit checks, and the classical fallback remain -part of the PhAST workflow. - -```yaml -solver: - damage_update: learned_replacement - damage_predictor: examples.learned_damage.predictor_plugin:create_predictor - damage_checkpoint: path/to/model.pt - damage_residual_rtol: 1.0e-3 - damage_residual_atol: 1.0e-8 - damage_bound_tolerance: 1.0e-8 - damage_fallback: true -``` - -Replacement is accepted only after shape, finite-value, damage-bound, -irreversibility, phase-field Dirichlet, and projected-residual checks. If a -check fails, the default route is an exact classical fallback. Disabling -fallback turns rejection into an error and should be reserved for debugging. - -These safeguards are not evidence that a learned model is accurate or -generalizes. Scientific claims require independent comparisons against matched -finite-element trajectories and complete reporting of acceptance, rejection, -and fallback events. - -## 7. A defensible student workflow - -1. Run an existing classical example without changing its formulation. -2. Modify one component at a time: mesh, material, loading, or energy split. -3. Examine the route report and retained displacement, damage, and history. -4. Generate training data with the constitutive and discretization choices - intended at inference. -5. Introduce a learned model through `learned_proposal`. -6. Compare the corrected trajectory with an independent classical reference. -7. Study `learned_replacement` only after fallback reporting is understood. - -The plug-in is architecture-neutral. Students and developers are invited to -propose adapters, tutorials, and audited examples. If an installation step, -configuration field, or solver message is unclear, open a GitHub issue with -the configuration, environment report, and minimal reproduction. - -Continue with [Heterogeneous material fields](05_heterogeneous_material_fields.md) -for an element-ordered material-field teaching problem, or return to the -[example gallery](../example-gallery.md) for complete classical workflows. diff --git a/docs/tutorial/05_heterogeneous_material_fields.md b/docs/tutorial/05_heterogeneous_material_fields.md index e1aaef4e4..c282456f1 100644 --- a/docs/tutorial/05_heterogeneous_material_fields.md +++ b/docs/tutorial/05_heterogeneous_material_fields.md @@ -81,5 +81,5 @@ See the public [`examples/heterogeneous_fields/README.md`](https://github.com/CE for parameters, outputs, and adaptation guidance. Next, use the [mesh-resolution diagnostic](notebook_mesh_resolution.ipynb) to -interpret `h/l0`, or return to [Modular fracture problems](03_modular_fem_and_learned_damage.md) -before introducing a learned damage proposal. +interpret `h/l0`, or return to [Modular fracture problems](03_modular_fem.md) +to review how the problem components fit together. diff --git a/docs/tutorial/06_first_research_study.md b/docs/tutorial/06_first_research_study.md index 860aa93b1..765468d9e 100644 --- a/docs/tutorial/06_first_research_study.md +++ b/docs/tutorial/06_first_research_study.md @@ -101,7 +101,6 @@ choice, retained histories, and animation back to one public configuration. | Quasi-static fracture | `examples/quasistatic/notched_holed_plate/config.yaml` or the example identified as closest to the target problem | Requires convergence, mesh, parameter, and comparison evidence for a new claim. | | Explicit dynamic fracture | `examples/dynamic/B2_kalthoff_winkler/config.yaml` or another curated dynamic example | Requires stable time-step and output-cost assessment. | | Elementwise heterogeneous properties | `examples/heterogeneous_fields/parameters.yaml` with its script-contract runner | A bounded teaching route, not a general heterogeneous YAML adapter or validated microstructure study. | -| Learned damage proposal | The experimental interface described in [Modular FEM and learned damage](03_modular_fem_and_learned_damage.md) | Replaces only the damage proposal; no trained checkpoint is distributed and classical acceptance/fallback remains authoritative. | Do not construct a new configuration by combining options from unrelated examples unless the capability matrix and schema identify that combination as @@ -252,8 +251,7 @@ A suitable first milestone is: > mesh, solver, and evidence limitations. New constitutive laws, unsupported interface formulations, automatic image-to- -mesh pipelines, three-dimensional fracture, a new learned architecture, or a -new execution adapter are software-development projects. They should be scoped +mesh pipelines, three-dimensional fracture, or a new execution adapter are software-development projects. They should be scoped separately from the first modelling milestone. ## Requesting review or assistance diff --git a/docs/tutorial/index.md b/docs/tutorial/index.md index 2ed0d1c5b..ac87046f7 100644 --- a/docs/tutorial/index.md +++ b/docs/tutorial/index.md @@ -69,7 +69,7 @@ checkout, and prints the resolved commit before installation. | [Retained-results notebook](notebook_retained_results.ipynb) | 10 min | Inspect checked-in Miehe result evidence without rerunning the full calculation. | | [Python API](../user_guide/python_api.md) | 10-15 min | Author a model with `phast.Problem` and understand the fluent method map. | | [Visual glossary](02_visual_glossary.md) | 10 min | Read the picture-first guide to AT1/AT2, energy splits, and `l0`. | -| [Modular FEM and learned damage](03_modular_fem_and_learned_damage.md) | 20 min | Assemble geometry, material, boundary conditions, fracture choices, solver routes, and an audited learned-damage plug-in. | +| [Modular fracture problems](03_modular_fem.md) | 20 min | Assemble geometry, material, boundary conditions, fracture choices, and solver routes. | | [Heterogeneous material fields](05_heterogeneous_material_fields.md) | 15 min | Define element-ordered `E(x)` and `Gc(x)` arrays and solve a bounded AT2 damage teaching problem. | | [From tutorial to first research study](06_first_research_study.md) | 20 min | Convert a modelling question into a bounded PhAST study with explicit decisions, preflight checks, pilot execution, outputs, and validation evidence. | | [YAML workflow](../user_guide/yaml_workflow.md) | 10 min | Run a public declarative configuration and understand the standard result directory. | @@ -89,8 +89,8 @@ checkout, and prints the resolved commit before installation. 5. Inspect the completed run with `phast.load_result(...)`. 6. Build a small model with `phast.Problem`. 7. Read the [visual glossary](02_visual_glossary.md) if the terminology feels abstract. -8. Read [Modular FEM and learned damage](03_modular_fem_and_learned_damage.md) - before introducing a learned damage proposal. +8. Read [Modular fracture problems](03_modular_fem.md) + before changing one problem component at a time. 9. Run [Heterogeneous material fields](05_heterogeneous_material_fields.md) before adapting a segmented or multiphase material map. 10. Use [From tutorial to first research study](06_first_research_study.md) to diff --git a/docs/tutorial/notebook_layered_dcb.ipynb b/docs/tutorial/notebook_layered_dcb.ipynb index edb9821b9..8abeebd62 100644 --- a/docs/tutorial/notebook_layered_dcb.ipynb +++ b/docs/tutorial/notebook_layered_dcb.ipynb @@ -894,7 +894,7 @@ "| Result directory already exists | Restart and run all cells for a new study identifier. Do not overwrite another parameter choice. |\n", "| OneDrive reports conflicts | Finish writing before synchronisation elsewhere. HDF5 reduces file count but does not make concurrent writers safe. |\n", "\n", - "This example does not establish a general interface, plasticity, three-dimensional composite, dynamic multi-material, or learned-damage workflow. A neural predictor cannot be inserted into this particular adapter merely by changing its solver string. Check the [capability matrix](../user_guide/capability_matrix.md) before selecting another route.\n", + "This example does not establish a general interface, plasticity, three-dimensional composite, or dynamic multi-material workflow. Check the [capability matrix](../user_guide/capability_matrix.md) before selecting another route.\n", "\n", "The small-strain assumption, diffuse crack length relative to the circular region, elementwise boundary approximation, and confinement by the outer material affect interpretation. A convincing animation alone does not settle those modelling questions.\n", "\n", diff --git a/docs/user_guide/capability_matrix.md b/docs/user_guide/capability_matrix.md index ee1df5081..8c37e9bd6 100644 --- a/docs/user_guide/capability_matrix.md +++ b/docs/user_guide/capability_matrix.md @@ -42,8 +42,6 @@ see the [plasticity and interface beta workflow](../supported_workflows/plastici | Coupled brittle PF + cohesive elements | Beta | Staggered AT2 matrix damage can be combined with zero-thickness cohesive-interface delamination. Calibrated structural validation remains incomplete. | | PF-CZM | Beta | Wu PF-CZM is available through `pf_model: PFCZM` for nonlinear forward damage solves with tensile-strength-calibrated rational degradation and uniaxial tests. Structural crack-growth validation is incomplete, and PF-plasticity-cohesive coupling is not supported. | | Coupled PF + plasticity + cohesive interfaces | Unsupported | This combined formulation is not available through the documented solver interface. | -| Learned damage proposal | Experimental | A user-supplied predictor may provide a projected initial guess; the classical damage solve remains authoritative. No trained model is distributed. | -| Audited learned damage replacement | Experimental | A predicted damage state may replace one damage solve only after bound, irreversibility, phase-field boundary-condition, and projected-residual checks. Rejection uses the classical fallback by default. | | 3D fracture | Unsupported | The documented fracture element pathways are two-dimensional. | | P2 / Q8 / Q9 element primitives | Scaffold | Shape functions, quadrature, and single-element stiffness tests exist for higher-order families; global solver dispatch is not supported. | | Native Q4 isotropic mechanics + AT2 damage | Beta | Structured Q4 mesh helpers, native Q4 mesh input, 2x2-Gauss isotropic mechanics, SciPy/MUMPS sparse-direct stiffness assembly, and matrix-free Q4 AT2 damage are tested. Q4 PF-CZM, AT1, and plasticity are not supported. | diff --git a/docs/user_guide/learned_damage.md b/docs/user_guide/learned_damage.md deleted file mode 100644 index e21955fdf..000000000 --- a/docs/user_guide/learned_damage.md +++ /dev/null @@ -1,214 +0,0 @@ -# Learned damage predictor interface - -## Status and claim boundary - -| Route | Public status | Accepted state | -|---|---|---| -| `classical` | Supported default | Classical phase-field damage solve | -| `learned_proposal` | Experimental interface | Classical solve initialized by a projected prediction | -| `learned_replacement` | Experimental interface | Prediction only after admissibility and projected-residual audits | -| Bundled trained model | Not provided | Users supply and validate their own checkpoint | - -The interface supports research integration. It does not establish the -accuracy, generalization, or speed of any neural architecture. - -## Standard computational route - -PhAST retains one standard finite-element workflow. Geometry definition, mesh -construction, material assignment, initial conditions, loading, boundary -conditions, mechanical equilibrium or time integration, history-field -evaluation, and result generation are independent of the selected damage -update. Only the implementation used to obtain the next phase-field damage -state is selectable. - -```{mermaid} -flowchart LR - A["Geometry, mesh, material,
loads and boundary conditions"] - B["Mechanical FEM update"] - C["Evaluate strain energy
and history field H"] - D{"Damage-update route"} - E["Classical phase-field
damage solve"] - F["Learned damage proposal"] - G["Projection and exact
damage correction"] - H["Admissibility and
projected-residual audit"] - I["Accepted damage state"] - J["Classical fallback"] - K["Advance time or
load increment"] - - A --> B - B --> C - C --> D - D -->|"classical"| E - D -->|"learned_proposal"| F - F --> G - G --> E - D -->|"learned_replacement"| F - F --> H - H -->|"accepted"| I - H -->|"rejected"| J - J --> E - E --> I - I --> K - K --> B - - classDef input fill:#f5f5f4,stroke:#57534e,color:#1c1917,stroke-width:1.2px; - classDef fem fill:#e7e5e4,stroke:#44403c,color:#1c1917,stroke-width:1.2px; - classDef learned fill:#e2e8f0,stroke:#475569,color:#0f172a,stroke-width:1.2px; - classDef audit fill:#f1f5f9,stroke:#64748b,color:#0f172a,stroke-width:1.2px; - classDef decision fill:#fafaf9,stroke:#57534e,color:#1c1917,stroke-width:1.4px; - - class A input; - class B,C,E,I,J,K fem; - class F learned; - class G,H audit; - class D decision; -``` - -### Meaning of replacement - -The optional learned component replaces only the classical solution of the -damage subproblem for an individual update. It does not replace: - -- the finite-element mesh or interpolation; -- the mechanical displacement solve; -- constitutive evaluation; -- the tensile/compressive energy decomposition; -- history-field construction; -- material degradation in the mechanical equations; -- loading, boundary conditions, or time integration; -- output, provenance, or route reporting. - -The phrase *learned replacement* therefore means replacement of the damage -subproblem solve by a predicted nodal damage state. The classical damage solver -remains available as the default route and as the fallback authority. - -### Recommended standard route - -For ordinary simulations, use `damage_update: classical`. For initial research -with a trained model, use `damage_update: learned_proposal`. In this mode the -prediction is projected onto the admissible set and used only as an initial -guess; the classical phase-field solver still determines the accepted damage -state. - -Use `damage_update: learned_replacement` only for an explicitly identified -experimental study. A prediction is accepted only after the configured -physical and residual audits. Failure of any audit invokes the classical -damage solve when fallback is enabled. - -## Configuration - -```yaml -solver: - damage_update: learned_proposal - damage_predictor: my_package.damage_adapter:create_predictor - damage_checkpoint: checkpoints/model.pt - damage_predictor_options: - representation: damage_increment - damage_residual_rtol: 1.0e-3 - damage_residual_atol: 1.0e-8 - damage_bound_tolerance: 1.0e-8 - damage_fallback: true -``` - -The factory is loaded from `module:factory` and called as: - -```python -factory(checkpoint=checkpoint, device=device, options=options) -``` - -It returns an object with: - -```python -predictor.name -predictor.predict(context) -``` - -`predict(context)` returns either a nodal tensor, interpreted as the next -damage field, or: - -```python -DamagePrediction( - damage=prediction, - representation="damage_increment", - diagnostics={"checkpoint": "..."}, -) -``` - -## State available to a plug-in - -`DamageStepContext` exposes: - -- step, time, and load factor; -- node coordinates and element connectivity; -- displacement and velocity; -- element and nodal history fields; -- previous accepted damage; -- material and formulation metadata; -- execution device and dtype. - -`canonical_node_features()` provides a small teaching representation. -`graph_edge_index()` provides dependency-free directed mesh edges. Neither is -mandatory. Architecture adapters may construct graph edge features, -node-element incidence operators, temporal windows, DeepONet branch/trunk -inputs, regular-grid projections, reduced coordinates, or their own normalized -features. - -This separation is deliberate. PhAST owns finite-element state and acceptance; -the plug-in owns architecture-specific preprocessing, checkpoint -interpretation, and inference. - -## Plug-and-play boundary - -The solver-facing interface is plug-and-play: installing a predictor through a -factory and selecting it in YAML does not require modification of the PhAST -mechanics or damage-solver source code. A conforming adapter receives -`DamageStepContext` and returns `DamagePrediction`. - -An arbitrary checkpoint is not automatically portable. Neural architectures -can differ in feature ordering, normalization, temporal history, graph -construction, finite-element incidence operators, grid projection, output -representation, and checkpoint serialization. These differences must be -declared and implemented by the model adapter. - -The reference TorchScript factory demonstrates the contract but is not a -registry of validated models. A future YAML-level model manifest may standardize -known adapters, normalization records, output representations, and checkpoint -metadata. Until such a registry is implemented and tested, the accurate -description is **plug-in compatible**, rather than universally -checkpoint-compatible. - -## Acceptance and fallback - -For `learned_proposal`, PhAST: - -1. checks shape and finite values; -2. projects damage to `[0, 1]`; -3. enforces irreversibility and phase-field Dirichlet values; -4. supplies the result as the classical solver initial guess. - -For `learned_replacement`, PhAST additionally checks: - -1. unprojected bound and irreversibility violations; -2. phase-field Dirichlet consistency; -3. the projected phase-field residual relative to the previous state; -4. an absolute residual floor. - -Rejected replacements use the classical damage solver when -`damage_fallback: true`. Route reports include calls, proposals, accepted -replacements, failures, and fallbacks. - -## Portability - -The predictor receives the device selected by PhAST, but its factory remains -responsible for verifying that the checkpoint and operators support that -device. CPU execution is the portable baseline. CUDA, MPS, and optional sparse -backends depend on the installed PyTorch build and platform-specific packages. - -Automatic route selection is deterministic according to availability; it is -not a guarantee of the fastest route for every mesh or machine. The public CI -workflow is configured to check the core CPU installation and configuration -path on Linux, macOS, and Windows. Optional accelerator and sparse-direct paths -require their own environment-specific evidence. - -See [Modular fracture problems and learned damage updates](../tutorial/03_modular_fem_and_learned_damage.md) -for the student workflow and `examples/learned_damage/` for a minimal factory. diff --git a/examples/learned_damage/README.md b/examples/learned_damage/README.md deleted file mode 100644 index 63f9a464b..000000000 --- a/examples/learned_damage/README.md +++ /dev/null @@ -1,146 +0,0 @@ -# Learned damage plug-in example - -PhAST exposes a small predictor protocol for research on learned damage -updates. The finite-element solver remains responsible for mechanics, -history-field construction, constraints, and acceptance. - -The example provides two adapters: - -- `PersistencePredictor` returns the previously converged damage field. It is - only an interface demonstration and is not a trained surrogate. -- `TorchScriptDamagePredictor` loads a user-supplied TorchScript checkpoint. - The model receives canonical nodal features and element connectivity. - -Add the following fields to a runnable PhAST YAML configuration: - -```yaml -solver: - damage_update: learned_proposal - damage_predictor: examples.learned_damage.predictor_plugin:create_predictor - damage_checkpoint: path/to/model.pt - damage_predictor_options: - representation: damage_increment -``` - -`learned_proposal` uses the prediction only as an initial guess for the -classical damage solve. The converged finite-element update remains -authoritative. - -`learned_replacement` is an experimental option. Before a prediction can -replace the classical update, PhAST checks its shape, finite values, physical -bounds, irreversibility, prescribed phase-field values, and projected damage -residual. A rejected prediction falls back to the classical solve by default. - -The canonical nodal feature columns are: - -```text -x, y, history field, previous damage, displacement x, displacement y -``` - -The canonical features are a teaching convenience, not a mandatory model -schema. A plug-in receives the complete `DamageStepContext` and can call -`context.graph_edge_index()` for mesh-graph models, use coordinates as a -DeepONet trunk, retain a temporal state buffer, project fields to a regular -grid, or construct node-element incidence operators. Architecture-specific -normalization and checkpoint schemas remain inside the plug-in. - -Use `representation: damage` when the model predicts the next field, or -`representation: damage_increment` when it predicts an increment relative to -the previous accepted damage. Both pass through the same PhAST projection and -audit. - -This interface does not prescribe a neural architecture or training -procedure. A predictor intended for scientific use must be trained and -evaluated on data consistent with the mesh representation, material model, -loading protocol, nondimensionalization, and fracture formulation used at -inference time. - -Python callers can install a predictor directly: - -```python -from examples.learned_damage import create_predictor - -predictor = create_predictor(checkpoint="path/to/model.pt", device=solver.device) -solver.set_damage_predictor(predictor, mode="learned_proposal") -``` - -Start with `learned_proposal`. Treat direct replacement as a separate -experimental study requiring independent finite-element comparisons and -complete reporting of all fallback events. - -## Swapping architectures - -The finite-element workflow does not change when a learned predictor is used. -Geometry construction, meshing, named region identification, boundary -conditions, loading, material properties, solver selection and time integration -are written exactly as for a classical run. One key inside the solver block -selects which operation performs the damage update, and a second names the -adapter: - -```yaml -solver: - solver_type: explicit - damage_update: learned_proposal - damage_predictor: examples.learned_damage.architectures.mesh_graph_net:create_predictor - damage_checkpoint: checkpoints/mesh_graph_net.pt - damage_fallback: true -``` - -Changing architecture means changing those two lines: - -| Architecture | `damage_predictor` | -|---|---| -| Persistence, interface demonstration | `examples.learned_damage.predictor_plugin:create_predictor` | -| TorchScript, architecture-agnostic | `examples.learned_damage.predictor_plugin:create_predictor` | -| Mesh-graph network | `examples.learned_damage.architectures.mesh_graph_net:create_predictor` | -| Your own | `my_package.my_architecture:create_predictor` | - -Each architecture needs a wrapper, because each one has its own graph -construction, feature ordering and normalization. The wrapper is the only place -that knowledge belongs. `examples/learned_damage/architectures/template.py` is a -commented skeleton with the two methods to fill in. - -## Finite-element helpers available to any predictor - -`DamageStepContext` supplies the quantities that predictors of every -architecture tend to need, so that a wrapper carries model-specific code only: - -| Method | Returns | -|---|---| -| `canonical_node_features()` | `[x, y, H, d_prev, u_x, u_y]` per node | -| `graph_edge_index()` | directed, duplicate-free mesh edges | -| `element_areas()` | area of every triangle | -| `assemble_element_to_nodes(v)` | `sum_e int N_i v_e dOmega` | -| `boundary_node_mask()` | topological boundary nodes, interior voids included | -| `edge_lengths(edge_index)` | length of every edge | - -## Portable checkpoints - -`architectures/mesh_graph_net.py` provides `export_torchscript`, which converts a -state-dict checkpoint into a TorchScript archive. The archive carries its own -architecture, so it evaluates with no model source on the import path and can be -distributed as a single file. Optimizer state is dropped in the process. - -```python -from examples.learned_damage.architectures.mesh_graph_net import export_torchscript - -export_torchscript("training_checkpoint.pt", "mesh_graph_net.pt") -``` - -Load checkpoints only from sources you trust. The state-dict loader accepts -tensor-and-metadata archives through PyTorch's restricted `weights_only` mode; -the example does not deserialize arbitrary Python objects. - -The `examples.learned_damage...` module paths are available when running from -the PhAST source checkout. They are examples rather than modules installed in -the core `phast` wheel. A deployed predictor should therefore live in its own -installed package or use an application-specific import path. - -## Auditing an accepted damage field - -`phast.solvers.damage_solver.damage_kkt_metrics` reports the box-constrained -KKT conditions for one damage field: projected stationarity, feasibility of the -irreversibility interval, dual sign conditions on the active sets, and -complementarity. Post-clamping and relaxed iterates can hide exactly these -violations, so the projected norms are the right quantity to report when -comparing a learned update against a classical one. diff --git a/examples/learned_damage/__init__.py b/examples/learned_damage/__init__.py deleted file mode 100644 index 0fee99b33..000000000 --- a/examples/learned_damage/__init__.py +++ /dev/null @@ -1,13 +0,0 @@ -"""Reference learned-damage plug-ins for PhAST.""" - -from .predictor_plugin import ( - PersistencePredictor, - TorchScriptDamagePredictor, - create_predictor, -) - -__all__ = [ - "PersistencePredictor", - "TorchScriptDamagePredictor", - "create_predictor", -] diff --git a/examples/learned_damage/architectures/__init__.py b/examples/learned_damage/architectures/__init__.py deleted file mode 100644 index 3fa627eb7..000000000 --- a/examples/learned_damage/architectures/__init__.py +++ /dev/null @@ -1,8 +0,0 @@ -"""Per-architecture learned-damage adapters. - -Each module here wraps one neural architecture so that it satisfies the PhAST -damage-predictor protocol. The solver, the finite-element workflow and the YAML -schema are identical for all of them; only the ``damage_predictor`` entry in the -solver block changes. Use :mod:`template` as the starting point for a new -architecture. -""" diff --git a/examples/learned_damage/architectures/mesh_graph_net.py b/examples/learned_damage/architectures/mesh_graph_net.py deleted file mode 100644 index 38d4eee84..000000000 --- a/examples/learned_damage/architectures/mesh_graph_net.py +++ /dev/null @@ -1,424 +0,0 @@ -"""Mesh-graph network adapter for the learned damage subproblem. - -One architecture, wrapped so that PhAST can call it. The finite-element -workflow is untouched: geometry, mesh, named regions, boundary conditions, -loading, material and time integration all remain exactly as they are for a -classical run. The only change is the damage subproblem's ``damage_predictor`` -entry. - -The generic finite-element quantities this adapter needs are supplied by -:class:`phast.learned_damage.DamageStepContext`. What is specific to this -architecture, and therefore lives here, is the feature contract: - -``node_features[:, 0]`` - ``sum_e int N_i H_e / (Gc * l0) dOmega``, the assembled and normalized - history integral. The network applies ``log10(1 + x)`` to this column - internally, so the raw assembled value is passed. -``node_features[:, 1]`` - A boundary indicator. See :func:`bounding_box_indicator` for the variant - used when the reference checkpoints were trained. -``edge_features[:, 0]`` - Edge length divided by ``l0``. - -A checkpoint evaluated under a different normalization, assembly or edge -feature still produces a plausible crack pattern while no longer representing -the trained operator. Reproducing the training-time construction exactly is the -adapter's responsibility, not the solver's. - -Usage ------ -.. code-block:: yaml - - solver: - damage_update: learned_proposal - damage_predictor: examples.learned_damage.architectures.mesh_graph_net:create_predictor - damage_checkpoint: checkpoints/mesh_graph_net.pt - damage_predictor_options: - representation: damage - boundary_indicator: bounding_box - damage_fallback: true -""" - -from __future__ import annotations - -from pathlib import Path -from typing import Any, Mapping - -import torch -from torch import nn - -from phast.learned_damage import DamagePrediction, DamageStepContext - -__all__ = [ - "MeshGraphDamageNet", - "MeshGraphNetPredictor", - "bounding_box_indicator", - "create_predictor", - "export_torchscript", - "load_network", -] - - -# ---------------------------------------------------------------- network - -def _mlp(in_features: int, out_features: int, hidden: int, layers: int) -> nn.Sequential: - modules: list[nn.Module] = [] - width = in_features - for _ in range(layers): - modules.append(nn.Linear(width, hidden)) - modules.append(nn.ReLU()) - width = hidden - modules.append(nn.Linear(width, out_features)) - return nn.Sequential(*modules) - - -class SymmetricProcessorBlock(nn.Module): - """One permutation-symmetric message-passing block.""" - - def __init__(self, latent_dim: int, hidden_layers: int) -> None: - super().__init__() - triple = 3 * latent_dim - self.edge_update = _mlp(triple, latent_dim, latent_dim, hidden_layers) - self.message_norm = nn.LayerNorm(triple) - self.message = _mlp(triple, latent_dim, latent_dim, hidden_layers) - self.node_norm = nn.LayerNorm(2 * latent_dim) - self.node_update = _mlp(2 * latent_dim, latent_dim, latent_dim, hidden_layers) - - def forward( - self, - node_latent: torch.Tensor, - edge_latent: torch.Tensor, - edge_index: torch.Tensor, - ) -> tuple[torch.Tensor, torch.Tensor]: - # Indexed rather than unpacked so that the block stays scriptable. - src, dst = edge_index[0], edge_index[1] - symmetric = torch.cat( - [ - node_latent[src] + node_latent[dst], - torch.abs(node_latent[src] - node_latent[dst]), - edge_latent, - ], - dim=-1, - ) - edge_latent = edge_latent + self.edge_update(symmetric) - symmetric = torch.cat( - [ - node_latent[src] + node_latent[dst], - torch.abs(node_latent[src] - node_latent[dst]), - edge_latent, - ], - dim=-1, - ) - messages = self.message(self.message_norm(symmetric)) - aggregate = torch.zeros_like(node_latent) - aggregate.index_add_(0, dst, messages) - degree = torch.zeros( - (node_latent.shape[0], 1), - dtype=node_latent.dtype, - device=node_latent.device, - ) - degree.index_add_( - 0, - dst, - torch.ones((dst.numel(), 1), dtype=node_latent.dtype, device=dst.device), - ) - aggregate = aggregate / degree.clamp_min(1.0) - node_input = self.node_norm(torch.cat([node_latent, aggregate], dim=-1)) - node_latent = node_latent + self.node_update(node_input) - return node_latent, edge_latent - - -class MeshGraphDamageNet(nn.Module): - """Encoder-processor-decoder mesh-graph surrogate for the damage field. - - The decoder is linear. Box constraints and irreversibility are applied by - the PhAST runtime coupling rather than by an output activation here. - """ - - def __init__( - self, - *, - latent_dim: int = 32, - hidden_layers: int = 3, - message_passing_steps: int = 10, - ) -> None: - super().__init__() - self.node_encoder = nn.Sequential( - _mlp(2, latent_dim, latent_dim, hidden_layers), - nn.LayerNorm(latent_dim), - ) - self.edge_encoder = nn.Sequential( - _mlp(1, latent_dim, latent_dim, hidden_layers), - nn.LayerNorm(latent_dim), - ) - self.processor = nn.ModuleList( - SymmetricProcessorBlock(latent_dim, hidden_layers) - for _ in range(message_passing_steps) - ) - self.decoder = _mlp(latent_dim, 1, latent_dim, hidden_layers) - - def forward( - self, - node_features: torch.Tensor, - edge_index: torch.Tensor, - edge_features: torch.Tensor, - ) -> torch.Tensor: - transformed = node_features.clone() - transformed[:, 0] = torch.log10(1.0 + transformed[:, 0].clamp_min(0.0)) - node_latent = self.node_encoder(transformed) - edge_latent = self.edge_encoder(edge_features) - for block in self.processor: - node_latent, edge_latent = block(node_latent, edge_latent, edge_index) - return self.decoder(node_latent).squeeze(-1) - - -# ---------------------------------------------------------------- checkpoints - -_ARCHITECTURE_KEYS = ("latent_dim", "hidden_layers", "message_passing_steps") - - -def load_network( - checkpoint: str | Path, - device: torch.device | str = "cpu", -) -> tuple[nn.Module, dict[str, Any]]: - """Load a TorchScript archive or a state-dict checkpoint. - - A TorchScript archive carries its own architecture, so evaluating it needs - no model source. A state-dict checkpoint must supply the three architecture - entries in a ``configuration`` mapping. - """ - path = Path(checkpoint) - try: - model = torch.jit.load(str(path), map_location=device) - model.eval() - return model, {"format": "torchscript", "path": str(path.resolve())} - except (RuntimeError, ValueError): - pass - - payload = torch.load(str(path), map_location="cpu", weights_only=True) - if not isinstance(payload, Mapping) or "model_state_dict" not in payload: - raise ValueError( - f"{path} is neither a TorchScript archive nor a checkpoint holding " - "'model_state_dict'") - configuration = dict(payload.get("configuration", {})) - missing = [key for key in _ARCHITECTURE_KEYS if key not in configuration] - if missing: - raise ValueError( - f"{path} is missing architecture entries {missing} in its " - "'configuration' mapping") - model = MeshGraphDamageNet( - latent_dim=int(configuration["latent_dim"]), - hidden_layers=int(configuration["hidden_layers"]), - message_passing_steps=int(configuration["message_passing_steps"]), - ) - model.load_state_dict(payload["model_state_dict"], strict=True) - model.to(device).eval() - metadata = { - "format": "state_dict", - "path": str(path.resolve()), - **{key: configuration[key] for key in _ARCHITECTURE_KEYS}, - } - for key in ("variant", "Gc", "length_scale", "seed"): - if key in configuration: - metadata[key] = configuration[key] - if "completed_epochs" in payload: - metadata["completed_epochs"] = int(payload["completed_epochs"]) - return model, metadata - - -def export_torchscript( - checkpoint: str | Path, - output: str | Path, - *, - n_nodes: int = 64, - verify: bool = True, -) -> Path: - """Convert a state-dict checkpoint into a portable TorchScript archive. - - The archive evaluates without this module on the import path, so trained - weights can be distributed without model source. Optimizer state is - dropped, which makes the archive substantially smaller than a training - checkpoint. - """ - model, metadata = load_network(checkpoint, device="cpu") - if metadata["format"] == "torchscript": - raise ValueError(f"{checkpoint} is already a TorchScript archive") - - scripted = torch.jit.script(model) - if verify: - generator = torch.Generator().manual_seed(0) - node_features = torch.rand(n_nodes, 2, generator=generator) - n_edges = 4 * n_nodes - edge_index = torch.randint(0, n_nodes, (2, n_edges), generator=generator) - edge_features = torch.rand(n_edges, 1, generator=generator) - with torch.no_grad(): - reference = model(node_features, edge_index, edge_features) - candidate = scripted(node_features, edge_index, edge_features) - assert torch.equal(reference, candidate), ( - "scripted model disagrees with the eager model") - - output = Path(output) - output.parent.mkdir(parents=True, exist_ok=True) - torch.jit.save(scripted, str(output)) - return output - - -# ---------------------------------------------------------------- features - -def bounding_box_indicator( - nodes: torch.Tensor, *, atol: float = 1.0e-10 -) -> torch.Tensor: - """Flag nodes on the bounding box of the mesh. - - This marks the outer extent only, so interior voids such as holes and slots - carry no flag. It is retained because the reference checkpoints were trained - against it. For the topological boundary, which does detect interior voids, - use :meth:`phast.learned_damage.DamageStepContext.boundary_node_mask`, - selected here by ``boundary_indicator: topological``. - """ - x, y = nodes[:, 0], nodes[:, 1] - on_edge = ( - torch.isclose(x, x.min(), atol=atol) - | torch.isclose(x, x.max(), atol=atol, rtol=0.0) - | torch.isclose(y, y.min(), atol=atol, rtol=0.0) - | torch.isclose(y, y.max(), atol=atol, rtol=0.0) - ) - return on_edge.to(dtype=nodes.dtype) - - -# ---------------------------------------------------------------- adapter - -class MeshGraphNetPredictor: - """Wrap :class:`MeshGraphDamageNet` for the PhAST damage-predictor protocol. - - Parameters - ---------- - model - A network taking ``(node_features, edge_index, edge_features)``. - representation - ``'damage'`` when the network predicts the next field, or - ``'damage_increment'`` when it predicts an increment relative to the - previous accepted damage. - boundary_indicator - ``'bounding_box'`` to reproduce the training-time construction, or - ``'topological'`` to mark interior voids as well. - length_scale, Gc - Override the values taken from the running material. Supply these only - when the checkpoint was trained under a different normalization, and - record the choice in any reported result. - metadata - Provenance recovered from the checkpoint, reported in diagnostics. - """ - - name = "mesh-graph-net" - - def __init__( - self, - model: nn.Module, - *, - representation: str = "damage", - boundary_indicator: str = "bounding_box", - length_scale: float | None = None, - Gc: float | None = None, - metadata: Mapping[str, Any] | None = None, - ) -> None: - if boundary_indicator not in {"bounding_box", "topological"}: - raise ValueError( - "boundary_indicator must be 'bounding_box' or 'topological'; " - f"got {boundary_indicator!r}") - self.model = model.eval() - self.representation = str(representation) - self.boundary_indicator = boundary_indicator - self.length_scale = length_scale - self.Gc = Gc - self.metadata = dict(metadata or {}) - self._graph_cache: tuple[ - tuple[int, int, int, int, float], torch.Tensor, torch.Tensor - ] | None = None - - @staticmethod - def _resolve(name: str, override: float | None, context: DamageStepContext) -> float: - if override is not None: - return float(override) - value = context.material.get(name) - if value is None: - raise ValueError( - f"the running material exposes no {name!r}; pass it explicitly " - "to MeshGraphNetPredictor") - return float(value) - - def _graph(self, context: DamageStepContext) -> tuple[torch.Tensor, torch.Tensor]: - """Edges and edge features, rebuilt only when the mesh changes.""" - length_scale = self._resolve("l0", self.length_scale, context) - key = ( - int(context.nodes.data_ptr()), - int(context.elements.data_ptr()), - int(context.nodes.shape[0]), - int(context.elements.shape[0]), - float(length_scale), - ) - if self._graph_cache is not None and self._graph_cache[0] == key: - return self._graph_cache[1], self._graph_cache[2] - edge_index = context.graph_edge_index() - edge_features = ( - context.edge_lengths(edge_index) / length_scale).reshape(-1, 1) - self._graph_cache = (key, edge_index, edge_features) - return edge_index, edge_features - - def node_features(self, context: DamageStepContext) -> torch.Tensor: - """Build the two nodal features this network was trained on.""" - Gc = self._resolve("Gc", self.Gc, context) - length_scale = self._resolve("l0", self.length_scale, context) - assembled = context.assemble_element_to_nodes( - context.history_element.reshape(-1) / (Gc * length_scale)) - if self.boundary_indicator == "bounding_box": - indicator = bounding_box_indicator(context.nodes) - else: - indicator = context.boundary_node_mask().to(dtype=assembled.dtype) - return torch.stack((assembled, indicator.to(assembled.dtype)), dim=1) - - def predict(self, context: DamageStepContext) -> DamagePrediction: - edge_index, edge_features = self._graph(context) - features = self.node_features(context) - model_dtype = next( - (parameter.dtype for parameter in self.model.parameters()), - torch.float32, - ) - with torch.no_grad(): - raw = self.model( - features.to(device=context.device, dtype=model_dtype), - edge_index.to(device=context.device), - edge_features.to(device=context.device, dtype=model_dtype), - ) - damage = torch.as_tensor( - raw, device=context.device, dtype=context.dtype, - ).reshape_as(context.damage_previous) - return DamagePrediction( - damage=damage, - representation=self.representation, - diagnostics={ - "adapter": self.name, - "n_edges": int(edge_index.shape[1]), - "boundary_indicator": self.boundary_indicator, - **self.metadata, - }, - ) - - -def create_predictor( - *, - checkpoint: str | Path, - device: torch.device | str = "cpu", - options: Mapping[str, Any] | None = None, -) -> MeshGraphNetPredictor: - """Factory matching the ``module:factory`` entry in a YAML configuration.""" - options = dict(options or {}) - model, metadata = load_network(checkpoint, device=device) - return MeshGraphNetPredictor( - model, - representation=options.get("representation", "damage"), - boundary_indicator=options.get("boundary_indicator", "bounding_box"), - length_scale=options.get("length_scale"), - Gc=options.get("Gc"), - metadata=metadata, - ) diff --git a/examples/learned_damage/architectures/template.py b/examples/learned_damage/architectures/template.py deleted file mode 100644 index f8e1b093f..000000000 --- a/examples/learned_damage/architectures/template.py +++ /dev/null @@ -1,141 +0,0 @@ -"""Template for a new learned-damage architecture. - -Copy this file, rename it after the architecture, and fill in the two marked -places. Nothing outside this file changes: the geometry, mesh, named regions, -boundary conditions, loading, material, solver and time integration are the same -as for a classical run, and the simulation selects the architecture through one -YAML entry. - -.. code-block:: yaml - - solver: - damage_update: learned_proposal - damage_predictor: my_package.my_architecture:create_predictor - damage_checkpoint: checkpoints/my_model.pt - damage_predictor_options: - representation: damage - damage_fallback: true - -What PhAST guarantees ---------------------- -A predictor is called once per damage subproblem and receives a -:class:`phast.learned_damage.DamageStepContext` describing the current -finite-element state. PhAST owns the admissibility audit, the projected -residual check, the box and irreversibility constraints, route reporting and -the classical fallback. A predictor proposes; it never decides. - -What the context offers ------------------------ -Raw state - ``nodes``, ``elements``, ``displacement``, ``velocity``, - ``history_element``, ``history_nodal``, ``damage_previous``, ``material``, - ``step``, ``time``, ``load_factor``, ``device``, ``dtype``. - -Finite-element helpers, architecture-neutral - ``canonical_node_features()`` returns ``[x, y, H, d_prev, u_x, u_y]``. - ``graph_edge_index()`` returns directed, duplicate-free mesh edges. - ``element_areas()`` returns the area of every triangle. - ``assemble_element_to_nodes(values)`` returns - ``sum_e int N_i value_e dOmega``. - ``boundary_node_mask()`` flags topological boundary nodes, interior voids - included. - ``edge_lengths(edge_index)`` returns the length of every edge. - -Anything beyond these, such as a radius graph, a regular-grid projection, a -trunk basis for an operator network, a temporal state buffer, or a -normalization recovered from a training manifest, is built here. That is the -whole purpose of a per-architecture wrapper: the solver stays neutral and each -model brings its own preprocessing. - -Obligations ------------ -1. Return one value per mesh node, shaped like ``context.damage_previous``. -2. Return it on ``context.device`` in ``context.dtype``. -3. Declare ``representation``: ``'damage'`` for the next field, or - ``'damage_increment'`` for an increment on the previous accepted damage. -4. Reproduce the training-time feature construction exactly. A mismatch in - normalization or feature order still yields a plausible crack pattern while - no longer representing the trained operator. -""" - -from __future__ import annotations - -from pathlib import Path -from typing import Any, Mapping - -import torch - -from phast.learned_damage import DamagePrediction, DamageStepContext - -__all__ = ["TemplatePredictor", "create_predictor"] - - -class TemplatePredictor: - """Minimal predictor skeleton. Replace the two marked bodies.""" - - #: Reported in the route summary and in prediction diagnostics. - name = "template" - - def __init__( - self, - model: Any, - *, - representation: str = "damage", - metadata: Mapping[str, Any] | None = None, - ) -> None: - self.model = model - self.representation = str(representation) - self.metadata = dict(metadata or {}) - - def build_inputs(self, context: DamageStepContext) -> tuple[Any, ...]: - """Build whatever this architecture consumes. - - Replace this body. The example below is the simplest useful case: the - canonical nodal features and the mesh graph. An operator network might - return coordinates and a history field instead; a radius-graph model - would build its own neighbourhood here. - """ - features = context.canonical_node_features() - edge_index = context.graph_edge_index() - return features, edge_index - - def predict(self, context: DamageStepContext) -> DamagePrediction: - """Call the model and return one nodal damage field. - - Replace the model call. Keep the surrounding contract: no gradient, one - value per node, the context's device and dtype. - """ - inputs = self.build_inputs(context) - with torch.no_grad(): - raw = self.model(*inputs) # replace with this architecture's call - - damage = torch.as_tensor( - raw, device=context.device, dtype=context.dtype, - ).reshape_as(context.damage_previous) - return DamagePrediction( - damage=damage, - representation=self.representation, - diagnostics={"adapter": self.name, **self.metadata}, - ) - - -def create_predictor( - *, - checkpoint: str | Path, - device: torch.device | str = "cpu", - options: Mapping[str, Any] | None = None, -) -> TemplatePredictor: - """Factory matching the ``module:factory`` entry in a YAML configuration. - - PhAST calls this with the checkpoint path, the resolved device and the - ``damage_predictor_options`` mapping. Loading, architecture reconstruction - and checkpoint validation belong here. - """ - options = dict(options or {}) - model = torch.jit.load(str(checkpoint), map_location=device) - model.eval() - return TemplatePredictor( - model, - representation=options.get("representation", "damage"), - metadata={"path": str(Path(checkpoint).resolve())}, - ) diff --git a/examples/learned_damage/predictor_plugin.py b/examples/learned_damage/predictor_plugin.py deleted file mode 100644 index b0508a0bf..000000000 --- a/examples/learned_damage/predictor_plugin.py +++ /dev/null @@ -1,75 +0,0 @@ -"""Minimal learned-damage plug-in contract. - -The persistence predictor is an interface demonstration, not a useful damage -model. The TorchScript adapter shows one portable way to load a trained model -without coupling PhAST to a particular neural-network architecture. -""" - -from __future__ import annotations - -from pathlib import Path -from typing import Any, Mapping, Optional - -import torch - -from phast.learned_damage import DamagePrediction, DamageStepContext - - -class PersistencePredictor: - """Return the previous converged damage field.""" - - name = "persistence-demonstration" - - def predict(self, context: DamageStepContext) -> DamagePrediction: - return DamagePrediction( - damage=context.damage_previous.clone(), - diagnostics={"purpose": "interface demonstration"}, - ) - - -class TorchScriptDamagePredictor: - """Adapt a TorchScript model to the PhAST damage-predictor protocol.""" - - def __init__( - self, - checkpoint: str, - device: torch.device, - representation: str = "damage", - ) -> None: - self.name = f"torchscript:{Path(checkpoint).name}" - self.representation = representation - self.model = torch.jit.load(checkpoint, map_location=device) - self.model.eval() - - def predict(self, context: DamageStepContext) -> DamagePrediction: - with torch.no_grad(): - damage = self.model( - context.canonical_node_features(), - context.elements, - ) - return DamagePrediction( - damage=torch.as_tensor( - damage, - device=context.device, - dtype=context.dtype, - ).reshape_as(context.damage_previous), - representation=self.representation, - diagnostics={"adapter": "torchscript"}, - ) - - -def create_predictor( - checkpoint: Optional[str] = None, - device: Optional[torch.device] = None, - options: Optional[Mapping[str, Any]] = None, -): - """Factory used by ``solver.damage_predictor`` configuration.""" - - options = dict(options or {}) - if checkpoint: - return TorchScriptDamagePredictor( - checkpoint=checkpoint, - device=device or torch.device("cpu"), - representation=str(options.get("representation", "damage")), - ) - return PersistencePredictor() diff --git a/examples/standard_workflow/dynamic_sent/config.yaml b/examples/standard_workflow/dynamic_sent/config.yaml index 1c2f62401..4693df442 100644 --- a/examples/standard_workflow/dynamic_sent/config.yaml +++ b/examples/standard_workflow/dynamic_sent/config.yaml @@ -55,7 +55,7 @@ analysis_steps: controls: {protocol: simple, num_steps: 4, dt: 1.0e-8, ramp_type: smooth_step, t_ramp: 2.0e-8} # 8. Select existing Velocity-Verlet mechanics and the classical AT2 damage solve. -# No quasi-static Newton backend, learned model, or three-dimensional route is used. +# No quasi-static Newton backend or three-dimensional route is used. solver: type: explicit device: cpu diff --git a/examples/standard_workflow/quasistatic_sent/config.yaml b/examples/standard_workflow/quasistatic_sent/config.yaml index aba582444..15464c6e8 100644 --- a/examples/standard_workflow/quasistatic_sent/config.yaml +++ b/examples/standard_workflow/quasistatic_sent/config.yaml @@ -54,7 +54,7 @@ analysis_steps: controls: {protocol: simple, num_steps: 4, dt: 0.25} # 8. Select existing staggered equilibrium/damage solves with SciPy mechanics. -# Tolerances control the existing solvers; no dynamics or learned model is used. +# Tolerances control the existing solvers; no dynamics is used. solver: type: quasi_static device: cpu diff --git a/src/phast/__init__.py b/src/phast/__init__.py index 0e20eaf02..fc7c1156e 100644 --- a/src/phast/__init__.py +++ b/src/phast/__init__.py @@ -67,11 +67,6 @@ DirectSolver, LBFGSSolver, MonolithicSolver, ) from .staggered_solver import StaggeredSolver, SolverConfig - from .learned_damage import ( - DamageDecision, DamagePrediction, DamagePredictionRejected, - DamagePredictor, DamageStepContext, DamageUpdateController, - load_damage_predictor, - ) from .device import DeviceContext, detect_device, get_device_tier, estimate_vram_mb, Profiler from .problem import Problem from .result import Result, ResultLoadError, load_result @@ -132,9 +127,6 @@ 'ExplicitDynamics', 'StaticSolver', 'QuasiStaticSolver', 'SecantCGSolver', 'LBFGSSolver', 'MonolithicSolver', 'StaggeredSolver', 'SolverConfig', - 'DamageDecision', 'DamagePrediction', 'DamagePredictionRejected', - 'DamagePredictor', 'DamageStepContext', 'DamageUpdateController', - 'load_damage_predictor', # Device / performance 'DeviceContext', 'detect_device', 'get_device_tier', 'estimate_vram_mb', 'Profiler', # Convenience diff --git a/src/phast/config/config.py b/src/phast/config/config.py index 2ae10210d..dc70387f6 100644 --- a/src/phast/config/config.py +++ b/src/phast/config/config.py @@ -379,16 +379,6 @@ class SolverSettings: softmax_H_beta: Optional[float] = None H_update_method: str = 'hard_max' enable_damage: bool = True - # Damage-subproblem route. ``classical`` is the default. Learned routes - # require a user-supplied ``module:factory`` predictor. - damage_update: str = 'classical' - damage_predictor: Optional[str] = None - damage_checkpoint: Optional[str] = None - damage_predictor_options: dict = field(default_factory=dict) - damage_residual_rtol: float = 1.0e-3 - damage_residual_atol: float = 1.0e-8 - damage_bound_tolerance: float = 1.0e-8 - damage_fallback: bool = True fail_on_mechanics_nonconvergence: bool = True fail_on_stagger_nonconvergence: bool = True adaptive_dt: bool = False @@ -1085,14 +1075,6 @@ def resolve_config(config: ProblemConfig) -> dict: softmax_H_beta=s.softmax_H_beta, H_update_method=s.H_update_method, enable_damage=s.enable_damage, - damage_update=s.damage_update, - damage_predictor=s.damage_predictor, - damage_checkpoint=s.damage_checkpoint, - damage_predictor_options=s.damage_predictor_options, - damage_residual_rtol=s.damage_residual_rtol, - damage_residual_atol=s.damage_residual_atol, - damage_bound_tolerance=s.damage_bound_tolerance, - damage_fallback=s.damage_fallback, fail_on_mechanics_nonconvergence=s.fail_on_mechanics_nonconvergence, fail_on_stagger_nonconvergence=s.fail_on_stagger_nonconvergence, adaptive_dt=s.adaptive_dt, diff --git a/src/phast/config/config_validation.py b/src/phast/config/config_validation.py index fed02b8f6..35eaa1b5c 100755 --- a/src/phast/config/config_validation.py +++ b/src/phast/config/config_validation.py @@ -78,9 +78,6 @@ # several anisotropic-PF papers since (issue #244). 'solver.stagger_norm': ['l2', 'linf'], 'solver.bounds_method': ['post_clamp', 'projected_cg'], - 'solver.damage_update': [ - 'classical', 'learned_proposal', 'learned_replacement', - ], # Issue #360 — H-update operator dispatcher. ``hard_max`` (default) # is byte-identical to ``torch.maximum``; the other methods are # opt-in differentiable alternatives. @@ -139,9 +136,6 @@ 'solver.H_cap_factor': (0.0, None), 'solver.stagger_tol': (0.0, None), 'solver.damage_tol': (0.0, None), - 'solver.damage_residual_rtol': (0.0, None), - 'solver.damage_residual_atol': (0.0, None), - 'solver.damage_bound_tolerance': (0.0, None), 'solver.static_tol': (0.0, None), 'material.sigma_ts': (0.0, None), 'material.pfczm_p': (2, None), @@ -800,8 +794,6 @@ def _validate_cross_field_compatibility(raw: dict, line_map: dict, solver_type = solver.get('solver_type', 'explicit') time_integrator = solver.get('time_integrator') - damage_update = solver.get('damage_update', 'classical') - damage_predictor = solver.get('damage_predictor') pf_model = str(_raw_material_value(raw, 'pf_model') or '').upper() degradation_type = str( _raw_material_value(raw, 'degradation_type') or 'standard' @@ -852,32 +844,6 @@ def _validate_cross_field_compatibility(raw: dict, line_map: dict, "damage-update route." ), )) - if damage_update != 'classical' and not damage_predictor: - errors.append(ValidationError( - path='solver.damage_predictor', - message=( - f"damage_update={damage_update!r} requires a predictor factory" - ), - line_no=line_map.get('solver.damage_update', - line_map.get('solver', 0)), - suggestion=( - "Set solver.damage_predictor to 'module:factory', or use " - "solver.damage_update: classical." - ), - )) - if damage_update == 'classical' and damage_predictor: - errors.append(ValidationError( - path='solver.damage_predictor', - message=( - "damage_predictor is configured but damage_update is classical" - ), - line_no=line_map.get('solver.damage_predictor', - line_map.get('solver', 0)), - suggestion=( - "Remove solver.damage_predictor or select " - "learned_proposal/learned_replacement explicitly." - ), - )) if solver_type != 'explicit' and time_integrator is not None: errors.append(ValidationError( path='solver.time_integrator', diff --git a/src/phast/core/fem_operators.py b/src/phast/core/fem_operators.py index aed655143..e0ebafdf5 100644 --- a/src/phast/core/fem_operators.py +++ b/src/phast/core/fem_operators.py @@ -1651,98 +1651,6 @@ def compute_crack_length(self, d, threshold=0.95): span = pts.max(0).values - pts.min(0).values return torch.norm(span).item() - def compute_physics_loss(self, u: torch.Tensor, d: torch.Tensor, - H_elem: torch.Tensor, bc_mask=None, - f_ext=None) -> torch.Tensor: - """Differentiable AT2 physics loss for neural operator training. - - Returns ||R_d||^2 + ||R_u||^2 where: - - R_d = AT2 damage residual (should be zero at equilibrium) - - R_u = mechanics residual f_int(u,d) - f_ext (should be zero at equilibrium) - - All operations preserve autograd graphs so gradients flow back - through u and d into the neural operator. - - Parameters - ---------- - u : (N, 2) displacement (from neural operator or solver) - d : (N,) damage (from neural operator or solver) - H_elem : (E,) driving force history variable - bc_mask : (N, 2) bool or None - If provided, constrained DOFs are excluded from the mechanics - residual (their equilibrium is enforced by BCs, not by R_u=0). - f_ext : (N, 2) or None - External force vector. If provided, the mechanics residual is - R_u = f_int - f_ext instead of R_u = f_int. - - Returns - ------- - loss : scalar tensor (differentiable) - """ - mesh = self.mesh - mat = self.material - pf_model = getattr(mat, 'pf_model', 'AT2').upper() - if pf_model != 'AT2': - raise NotImplementedError( - "FEMOperators.compute_physics_loss() currently implements " - f"the AT2 damage residual only; got pf_model={pf_model!r}. " - "Use AT2 or add a model-specific residual before using this " - "loss for AT1/PF-CZM training." - ) - Gc = mat.Gc - l0 = mat.l0 - - # --- Mechanics residual: f_int - f_ext should be zero at equilibrium --- - f_int = self.internal_force(u, d) - if f_ext is not None: - R_u = f_int - f_ext - else: - R_u = f_int - if bc_mask is not None: - R_u = R_u.clone() - R_u[bc_mask] = 0.0 # exclude constrained DOFs - R_u_sq = (R_u ** 2).sum() - - # --- Damage residual: R_d = A*d - b --- - # AT2 weak form: Gc*l0*K_lap*d + (2H + Gc/l0)*M*d - 2H*M*1 = 0 - # where M is the consistent mass matrix: area/12 * [2,1,1;1,2,1;1,1,2]. - # This matches PhaseFieldDamageSolver._Ax() exactly. - - # Laplacian stiffness term: Gc*l0 * K_lap @ d - K_lap_d = self.laplacian_matvec(d) - stiff_term = Gc * l0 * K_lap_d - - # Reaction (mass) term: (2H + Gc/l0) * M * d - # Consistent mass: M_e @ d_e = (area/12) * (d_a + d_sum) for node a - # where d_sum = sum of all 3 nodal d values in the element. - areas = mesh.areas - elements = mesh.elements - elem_flat = self._elem_flat - - rc = (2.0 * H_elem + Gc / l0) * areas / 12.0 # (E,) - d_e = d[elements] # (E, 3) - d_sum = d_e.sum(1) # (E,) - # Per-node contribution: rc * (d_a + d_sum) - react_contrib = rc.unsqueeze(1) * (d_e + d_sum.unsqueeze(1)) # (E, 3) - react_term = torch.zeros(mesh.n_nodes, dtype=self.dtype, - device=self.device) - react_term.scatter_add_(0, elem_flat, react_contrib.flatten()) - - # RHS: b = 2H * M * 1 (consistent mass applied to ones vector) - # M*1_a = area/12 * (1 + 3) = area/3 per node (since 1_sum = 3) - # But with nodal H: b_a = area/6 * (H_a + H_sum) — see damage_solver - # For element-level H, all 3 nodes get same H_e, so: - # b_a = 2*H_e * area/12 * (1 + 3) = 2*H_e * area/3 - rhs_coeff = 2.0 * H_elem * areas / 3.0 # (E,) - rhs_contrib = rhs_coeff.unsqueeze(1).expand(-1, 3).flatten() - b = torch.zeros(mesh.n_nodes, dtype=self.dtype, device=self.device) - b.scatter_add_(0, elem_flat, rhs_contrib) - - R_d = stiff_term + react_term - b - R_d_sq = (R_d ** 2).sum() - - return R_u_sq + R_d_sq - def compute_crack_front_indicator(self, d: torch.Tensor, d_prev: torch.Tensor, dt: float) -> torch.Tensor: """Compute crack front indicator as d * (dd/dt). diff --git a/src/phast/learned_damage.py b/src/phast/learned_damage.py deleted file mode 100644 index c1fff7947..000000000 --- a/src/phast/learned_damage.py +++ /dev/null @@ -1,550 +0,0 @@ -"""Public contracts for learned damage updates in staggered FEM simulations. - -The classical phase-field damage solver remains the default and authoritative -route. A learned model may be used in one of two explicitly selected modes: - -``learned_proposal`` - The prediction is projected to the admissible interval and supplied as the - initial guess to the classical damage solve. - -``learned_replacement`` - The prediction replaces one damage subproblem only when it passes finite - value, shape, bound, irreversibility, phase-field Dirichlet, and projected - residual checks. A rejected prediction falls back to the classical solve - unless fallback is disabled. - -This module defines inference contracts. It does not define a network -architecture, training procedure, dataset, or generalization claim. -""" - -from __future__ import annotations - -from dataclasses import dataclass, field -import importlib -import math -from typing import Any, Callable, Mapping, Protocol, runtime_checkable - -import torch - - -LEARNED_DAMAGE_MODES = { - "classical", - "learned_proposal", - "learned_replacement", -} - - -@dataclass(frozen=True) -class DamageStepContext: - """Read-only description of one damage-subproblem call. - - Tensor fields are detached views or clones supplied for inference. A model - must return one nodal damage value per mesh node. - """ - - step: int - time: float | None - load_factor: float - nodes: torch.Tensor - elements: torch.Tensor - displacement: torch.Tensor - velocity: torch.Tensor - history_element: torch.Tensor - history_nodal: torch.Tensor - damage_previous: torch.Tensor - material: Mapping[str, Any] - phase_field_model: str - energy_split: str - device: torch.device - dtype: torch.dtype - - def canonical_node_features(self) -> torch.Tensor: - """Return ``[x, y, H, d_prev, u_x, u_y]`` node features. - - The helper is intentionally small and architecture-neutral. Models that - require graph edges, quadrature history, normalization, or additional - material channels should construct them from the full context. - """ - - coordinates = self.nodes[:, :2].to( - device=self.device, dtype=self.dtype) - return torch.cat( - ( - coordinates, - self.history_nodal.reshape(-1, 1), - self.damage_previous.reshape(-1, 1), - self.displacement[:, :2], - ), - dim=1, - ) - - def graph_edge_index(self) -> torch.Tensor: - """Return directed, duplicate-free mesh edges as ``[2, n_edges]``. - - Plug-ins may ignore this helper, construct node-element incidence - operators, or add edge features. PhAST does not require PyTorch - Geometric for the public predictor contract. - """ - - local_nodes = self.elements.shape[1] - pairs = [] - for local_index in range(local_nodes): - left = self.elements[:, local_index] - right = self.elements[:, (local_index + 1) % local_nodes] - pairs.append(torch.stack((left, right), dim=1)) - pairs.append(torch.stack((right, left), dim=1)) - edges = torch.cat(pairs, dim=0) - return torch.unique(edges, dim=0).t().contiguous() - - def element_areas(self) -> torch.Tensor: - """Return the area of every three-node triangle in the mesh.""" - - if self.elements.shape[1] != 3: - raise NotImplementedError( - "element_areas currently covers three-node triangles; got " - f"{self.elements.shape[1]} local nodes") - xy = self.nodes[:, :2][self.elements] - x0, x1, x2 = xy[:, 0], xy[:, 1], xy[:, 2] - twice_area = ( - (x1[:, 0] - x0[:, 0]) * (x2[:, 1] - x0[:, 1]) - - (x2[:, 0] - x0[:, 0]) * (x1[:, 1] - x0[:, 1]) - ) - return 0.5 * twice_area.abs() - - def assemble_element_to_nodes(self, values: torch.Tensor) -> torch.Tensor: - """Assemble ``sum_e int N_i value_e dOmega`` over the mesh. - - This is the finite-element integral of a linear shape function against - a piecewise-constant element field, without mass inversion. Divide by - the lumped mass returned for a unit field to obtain a nodal average - instead. - """ - - values = values.reshape(-1) - if values.shape[0] != self.elements.shape[0]: - raise ValueError( - "values must hold one entry per element; got " - f"{values.shape[0]} for {self.elements.shape[0]} elements") - n_local = self.elements.shape[1] - weighted = ( - values * self.element_areas() / n_local - ).repeat_interleave(n_local) - accumulated = torch.zeros( - self.nodes.shape[0], dtype=weighted.dtype, device=weighted.device) - accumulated.index_add_(0, self.elements.reshape(-1), weighted) - return accumulated - - def boundary_node_mask(self) -> torch.Tensor: - """Flag nodes on the topological boundary of the mesh. - - An element edge shared by exactly one element lies on a boundary. The - test is topological, so interior voids such as holes and slots are - detected alongside the outer perimeter. - """ - - n_local = self.elements.shape[1] - pairs = torch.cat( - [ - torch.stack( - ( - self.elements[:, local], - self.elements[:, (local + 1) % n_local], - ), - dim=1, - ) - for local in range(n_local) - ], - dim=0, - ) - undirected, counts = torch.unique( - pairs.sort(dim=1).values, dim=0, return_counts=True) - boundary_edges = undirected[counts == 1] - mask = torch.zeros( - self.nodes.shape[0], dtype=torch.bool, device=self.nodes.device) - mask[boundary_edges.reshape(-1)] = True - return mask - - def edge_lengths(self, edge_index: torch.Tensor) -> torch.Tensor: - """Return the Euclidean length of every edge in ``edge_index``.""" - - delta = ( - self.nodes[:, :2][edge_index[0]] - self.nodes[:, :2][edge_index[1]] - ) - return torch.linalg.vector_norm(delta, dim=1) - - -@dataclass(frozen=True) -class DamagePrediction: - """Prediction returned by a learned damage model.""" - - damage: torch.Tensor - diagnostics: Mapping[str, Any] = field(default_factory=dict) - representation: str = "damage" - - -@runtime_checkable -class DamagePredictor(Protocol): - """Protocol implemented by a learned damage predictor.""" - - name: str - - def predict( - self, - context: DamageStepContext, - ) -> DamagePrediction | torch.Tensor: - """Return a nodal damage field or damage increment.""" - - -class DamagePredictionRejected(RuntimeError): - """Raised when replacement is rejected and fallback is disabled.""" - - -@dataclass(frozen=True) -class DamageDecision: - """Audited decision for one learned damage prediction.""" - - route: str - candidate: torch.Tensor | None - accepted_replacement: bool - reason: str - residual_relative: float | None = None - residual_rms: float | None = None - lower_bound_violation: float | None = None - upper_bound_violation: float | None = None - diagnostics: Mapping[str, Any] = field(default_factory=dict) - - -def _predictor_name(predictor: Any) -> str: - return str( - getattr( - predictor, - "name", - getattr(predictor, "__name__", predictor.__class__.__name__), - ) - ) - - -def load_damage_predictor( - specification: str, - *, - checkpoint: str | None, - device: torch.device, - options: Mapping[str, Any] | None = None, -) -> DamagePredictor | Callable[[DamageStepContext], Any]: - """Load a predictor from a public ``module:factory`` specification. - - The factory is called as:: - - factory(checkpoint=..., device=..., options={...}) - - It must return an object with ``predict(context)`` or a callable accepting - one :class:`DamageStepContext`. Checkpoint interpretation belongs to the - user-provided factory so PhAST does not impose one neural architecture. - """ - - if ":" not in specification: - raise ValueError( - "damage_predictor must use 'module:factory' notation; " - f"received {specification!r}" - ) - module_name, factory_name = specification.split(":", 1) - module = importlib.import_module(module_name) - factory = getattr(module, factory_name, None) - if factory is None or not callable(factory): - raise ValueError( - f"Damage predictor factory {factory_name!r} was not found as a " - f"callable in module {module_name!r}." - ) - predictor = factory( - checkpoint=checkpoint, - device=device, - options=dict(options or {}), - ) - if isinstance(predictor, torch.nn.Module): - predictor.eval() - if not hasattr(predictor, "predict") and not callable(predictor): - raise TypeError( - "A damage predictor factory must return an object with " - "predict(context) or a callable accepting DamageStepContext." - ) - return predictor - - -class DamageUpdateController: - """Audit learned predictions and select the damage-update route.""" - - def __init__( - self, - predictor: DamagePredictor | Callable[[DamageStepContext], Any], - *, - mode: str = "learned_proposal", - residual_rtol: float = 1.0e-3, - residual_atol: float = 1.0e-8, - bound_tolerance: float = 1.0e-8, - fallback: bool = True, - ) -> None: - if mode not in LEARNED_DAMAGE_MODES - {"classical"}: - raise ValueError( - "Learned damage mode must be 'learned_proposal' or " - f"'learned_replacement'; received {mode!r}." - ) - if residual_rtol < 0 or residual_atol < 0 or bound_tolerance < 0: - raise ValueError("Damage audit tolerances must be non-negative.") - self.predictor = predictor - self.mode = mode - self.residual_rtol = float(residual_rtol) - self.residual_atol = float(residual_atol) - self.bound_tolerance = float(bound_tolerance) - self.fallback = bool(fallback) - self.calls = 0 - self.proposals = 0 - self.accepted_replacements = 0 - self.fallbacks = 0 - self.failures = 0 - self.last_decision: DamageDecision | None = None - - @property - def predictor_name(self) -> str: - return _predictor_name(self.predictor) - - def _predict(self, context: DamageStepContext) -> DamagePrediction: - with torch.no_grad(): - if hasattr(self.predictor, "predict"): - output = self.predictor.predict(context) - else: - output = self.predictor(context) - if isinstance(output, DamagePrediction): - return output - if isinstance(output, torch.Tensor): - return DamagePrediction(output) - raise TypeError( - "Damage predictor output must be a torch.Tensor or " - "DamagePrediction." - ) - - @staticmethod - def _projected_residual( - residual: torch.Tensor, - damage: torch.Tensor, - damage_previous: torch.Tensor, - tolerance: float, - ) -> torch.Tensor: - lower_active = damage <= damage_previous + tolerance - upper_active = damage >= 1.0 - tolerance - projected = torch.where( - lower_active, - torch.minimum(residual, torch.zeros_like(residual)), - residual, - ) - projected = torch.where( - upper_active, - torch.maximum(projected, torch.zeros_like(projected)), - projected, - ) - return projected - - def _reject( - self, - reason: str, - *, - diagnostics: Mapping[str, Any] | None = None, - **metrics: Any, - ) -> DamageDecision: - self.failures += 1 - if self.fallback: - self.fallbacks += 1 - decision = DamageDecision( - route="classical_fallback", - candidate=None, - accepted_replacement=False, - reason=reason, - diagnostics=dict(diagnostics or {}), - **metrics, - ) - self.last_decision = decision - return decision - raise DamagePredictionRejected(reason) - - def decide( - self, - context: DamageStepContext, - *, - damage_solver: Any, - phase_field_mask: torch.Tensor | None = None, - phase_field_values: torch.Tensor | None = None, - ) -> DamageDecision: - """Return an audited proposal, replacement, or fallback decision.""" - - self.calls += 1 - try: - prediction = self._predict(context) - except Exception as exc: - return self._reject(f"predictor error: {exc}") - - diagnostics = dict(prediction.diagnostics) - diagnostics.setdefault("representation", prediction.representation) - raw_prediction = prediction.damage.detach().to( - device=context.device, dtype=context.dtype) - if raw_prediction.shape != context.damage_previous.shape: - return self._reject( - "prediction shape does not match nodal damage shape", - diagnostics=diagnostics, - ) - if prediction.representation == "damage": - candidate = raw_prediction - elif prediction.representation == "damage_increment": - candidate = context.damage_previous + raw_prediction - else: - return self._reject( - "prediction representation must be 'damage' or " - "'damage_increment'", - diagnostics=diagnostics, - ) - if not bool(torch.isfinite(candidate).all()): - return self._reject( - "prediction contains NaN or infinite values", - diagnostics=diagnostics, - ) - - damage_previous = context.damage_previous - lower_violation = float( - torch.clamp(damage_previous - candidate, min=0).max().item()) - upper_violation = float( - torch.clamp(candidate - 1.0, min=0).max().item()) - - if self.mode == "learned_proposal": - candidate = torch.maximum( - torch.clamp(candidate, min=0.0, max=1.0), - damage_previous, - ) - if phase_field_mask is not None and phase_field_values is not None: - candidate = torch.where( - phase_field_mask.to(device=context.device), - phase_field_values.to( - device=context.device, dtype=context.dtype), - candidate, - ) - self.proposals += 1 - decision = DamageDecision( - route="learned_proposal_exact_correction", - candidate=candidate, - accepted_replacement=False, - reason="admissible prediction supplied as exact-solver initial guess", - lower_bound_violation=lower_violation, - upper_bound_violation=upper_violation, - diagnostics=diagnostics, - ) - self.last_decision = decision - return decision - - if ( - lower_violation > self.bound_tolerance - or upper_violation > self.bound_tolerance - ): - return self._reject( - "prediction violates damage bounds or irreversibility", - diagnostics=diagnostics, - lower_bound_violation=lower_violation, - upper_bound_violation=upper_violation, - ) - - if phase_field_mask is not None and phase_field_values is not None: - mask = phase_field_mask.to(device=context.device) - values = phase_field_values.to( - device=context.device, dtype=context.dtype) - if bool(mask.any()): - bc_error = float( - torch.abs(candidate[mask] - values[mask]).max().item()) - if bc_error > self.bound_tolerance: - return self._reject( - "prediction violates a phase-field Dirichlet condition", - diagnostics=diagnostics, - lower_bound_violation=lower_violation, - upper_bound_violation=upper_violation, - ) - - try: - residual = damage_solver.compute_residual( - context.history_element, candidate) - reference = damage_solver.compute_residual( - context.history_element, damage_previous) - except Exception as exc: - return self._reject( - f"damage residual audit failed: {exc}", - diagnostics=diagnostics, - lower_bound_violation=lower_violation, - upper_bound_violation=upper_violation, - ) - - projected = self._projected_residual( - residual, - candidate, - damage_previous, - self.bound_tolerance, - ) - if phase_field_mask is not None: - projected = torch.where( - phase_field_mask.to(device=projected.device), - torch.zeros_like(projected), - projected, - ) - projected_norm = float(torch.linalg.vector_norm(projected).item()) - reference_norm = float(torch.linalg.vector_norm(reference).item()) - residual_relative = projected_norm / max( - reference_norm, - torch.finfo(projected.dtype).eps, - ) - residual_rms = projected_norm / math.sqrt(max(projected.numel(), 1)) - residual_ok = ( - residual_relative <= self.residual_rtol - or residual_rms <= self.residual_atol - ) - if not residual_ok: - return self._reject( - "prediction failed the projected damage-residual audit", - diagnostics=diagnostics, - residual_relative=residual_relative, - residual_rms=residual_rms, - lower_bound_violation=lower_violation, - upper_bound_violation=upper_violation, - ) - - self.accepted_replacements += 1 - decision = DamageDecision( - route="learned_replacement", - candidate=candidate.clone(), - accepted_replacement=True, - reason="prediction passed admissibility and residual audits", - residual_relative=residual_relative, - residual_rms=residual_rms, - lower_bound_violation=lower_violation, - upper_bound_violation=upper_violation, - diagnostics=diagnostics, - ) - self.last_decision = decision - return decision - - def summary(self) -> dict[str, Any]: - """Return machine-readable controller statistics.""" - - return { - "mode": self.mode, - "predictor": self.predictor_name, - "calls": self.calls, - "proposals": self.proposals, - "accepted_replacements": self.accepted_replacements, - "fallbacks": self.fallbacks, - "failures": self.failures, - "residual_rtol": self.residual_rtol, - "residual_atol": self.residual_atol, - "bound_tolerance": self.bound_tolerance, - "fallback_enabled": self.fallback, - "last_route": ( - None if self.last_decision is None - else self.last_decision.route - ), - "last_reason": ( - None if self.last_decision is None - else self.last_decision.reason - ), - } diff --git a/src/phast/solvers/damage_solver.py b/src/phast/solvers/damage_solver.py index 40e5ea723..31bdd30ec 100644 --- a/src/phast/solvers/damage_solver.py +++ b/src/phast/solvers/damage_solver.py @@ -2352,8 +2352,8 @@ def _solve_forward_with_material_overrides( pf_dirichlet_mask=None, pf_dirichlet_values=None): """Forward-only solve with temporary material overrides. - This path is used by no-grad dataset generation and live hybrid - audits. It temporarily installs the per-element material arrays, runs + This path is used by no-grad dataset generation and forward-only + residual checks. It temporarily installs the per-element material arrays, runs the normal constrained solver, and restores the cached scalar or gamma-corrected arrays afterwards. """ diff --git a/src/phast/solvers/staggered_solver.py b/src/phast/solvers/staggered_solver.py index d983d9d3a..a50be3262 100644 --- a/src/phast/solvers/staggered_solver.py +++ b/src/phast/solvers/staggered_solver.py @@ -16,7 +16,7 @@ - H capping for out-of-distribution stability (Pantidis 2026) - Modular sub-solvers (swappable mechanics + damage) - Support for quasi-static NR (not just explicit dynamics) - - Hooks for neural operator injection (set_damage, get_H_nodal) + - Hooks for external state access (set_damage, get_H_nodal) """ import torch @@ -39,11 +39,6 @@ init_h5, CSVHistory) from ..utils.device import DeviceContext from ..plasticity import DuctilePhaseFieldCoupling, MeshJ2Elastoplasticity -from ..learned_damage import ( - DamageStepContext, - DamageUpdateController, - load_damage_predictor, -) @dataclass @@ -112,14 +107,6 @@ class SolverConfig: # ExplicitDynamics docstring. Zero = pure Verlet (default). damping_ratio_max: float = 0.0 enable_damage: bool = True # False = pure elastic (skip damage solve + H update) - damage_update: str = 'classical' # classical | learned_proposal | learned_replacement - damage_predictor: Optional[str] = None # module:factory - damage_checkpoint: Optional[str] = None - damage_predictor_options: dict | None = None - damage_residual_rtol: float = 1.0e-3 - damage_residual_atol: float = 1.0e-8 - damage_bound_tolerance: float = 1.0e-8 - damage_fallback: bool = True # Linear-solver backend for the QuasiStaticSolver Newton inner step. # Forwarded to QuasiStaticSolver(backend=...) in _build_mechanics_solver. # 'auto' picks the best available sparse-direct path; 'scipy', 'mumps', @@ -168,18 +155,6 @@ class SolverConfig: print_every: int = 10 def __post_init__(self): - if self.damage_update not in { - 'classical', 'learned_proposal', 'learned_replacement', - }: - raise ValueError( - "damage_update must be classical, learned_proposal, or " - f"learned_replacement; got {self.damage_update!r}" - ) - if self.damage_update != 'classical' and not self.damage_predictor: - raise ValueError( - f"damage_update={self.damage_update!r} requires " - "damage_predictor='module:factory'." - ) if (self.preconditioner is None and self.solver_type in ('static', 'quasi_static', 'quasi_static_legacy', 'lbfgs', @@ -405,23 +380,6 @@ def __init__(self, mesh: FEMMesh, material: Material, self._last_mechanics_converged = True self._last_mechanics_iter = 0 self._last_dt_used = 0.0 # for adaptive-dt history (forward compat) - self.damage_controller = None - self._last_damage_route = 'classical' - if config.damage_update != 'classical': - predictor = load_damage_predictor( - config.damage_predictor, - checkpoint=config.damage_checkpoint, - device=self.device, - options=config.damage_predictor_options, - ) - self.set_damage_predictor( - predictor, - mode=config.damage_update, - residual_rtol=config.damage_residual_rtol, - residual_atol=config.damage_residual_atol, - bound_tolerance=config.damage_bound_tolerance, - fallback=config.damage_fallback, - ) print(f"[StaggeredSolver] State allocated: {N} nodes, {E} elements, " f"device={self.device}, dtype={self.dtype}", flush=True) self.print_route_report() @@ -925,15 +883,15 @@ def step_compute_driving_force(self, strain=None) -> torch.Tensor: self.H_elem = torch.clamp(self.H_elem, max=self.config.H_cap) # Note: elem_to_node performs area-weighted averaging, which smooths - # peak H values near the crack tip. Neural operators using get_H_nodal() - # will see a smoother driving force than the element-level H_elem that - # the damage solver actually uses. For NN training, consider using - # H_elem directly with element-to-node projection at inference time. + # peak H values near the crack tip. Callers of get_H_nodal() therefore + # see a smoother driving force than the element-level H_elem that the + # damage solver actually uses; read H_elem directly when the + # element-level values are required. self.H_nodal = self.mesh.elem_to_node(self.H_elem) return psi def step_solve_damage(self, d_prev_step=None): - """Update damage through the configured classical or learned route. + """Solve AT2 for damage. If `self.diff_Gc` and/or `self.diff_l0` are set (0-d tensors, typically with requires_grad=True), they are forwarded to the @@ -974,28 +932,10 @@ def step_solve_damage(self, d_prev_step=None): pf_mask, pf_vals = self.bcs.get_pf_dirichlet_mask_values() d_prev_arg = d_prev_step if d_prev_step is not None else self.d initial_guess = self.d if d_prev_step is not None else None - decision = None - if self.damage_controller is not None: - context = self._damage_step_context(d_prev_arg) - decision = self.damage_controller.decide( - context, - damage_solver=self.damage_solver, - phase_field_mask=pf_mask, - phase_field_values=pf_vals, - ) - self._last_damage_route = decision.route - if decision.accepted_replacement: - self.d = decision.candidate - self._apply_pf_dirichlet() - return - if decision.candidate is not None: - initial_guess = decision.candidate self.d = self.damage_solver.solve( self.H_elem, d_prev_arg, Gc=Gc, l0=l0, Gc_field=Gc_field, pf_dirichlet_mask=pf_mask, pf_dirichlet_values=pf_vals, initial_guess=initial_guess) - if decision is None: - self._last_damage_route = 'classical' # Phase-field Dirichlet enforcement (issue #213): re-pin the # damage value at every step on nodes flagged with the # pf_dirichlet BC (matches COMSOL pre-crack convention). @@ -1404,7 +1344,7 @@ def step_full(self) -> torch.Tensor: return psi # ------------------------------------------------------------------ # - # Neural operator integration hooks + # External state and material-field hooks # ------------------------------------------------------------------ # def install_diff_E_field(self, E_field: torch.Tensor): @@ -1467,86 +1407,21 @@ def install_diff_E_field(self, E_field: torch.Tensor): self._E_field_first_install = True def set_damage(self, d_new: torch.Tensor): - """Inject external damage field (e.g., from GINO).""" + """Replace the current damage field with an externally supplied one.""" self.d = d_new.clone() def get_H_nodal(self) -> torch.Tensor: - """Get current history variable at nodes (for neural operator input).""" + """Get the current history variable projected to the nodes.""" return self.H_nodal - def set_damage_predictor( - self, predictor, *, mode='learned_proposal', - residual_rtol=1.0e-3, residual_atol=1.0e-8, - bound_tolerance=1.0e-8, fallback=True): - """Install a learned damage predictor for subsequent damage updates. - - ``learned_proposal`` always retains the classical damage solve. - ``learned_replacement`` accepts a prediction only after the audits - implemented by :class:`phast.DamageUpdateController`. - """ - self.damage_controller = DamageUpdateController( - predictor, - mode=mode, - residual_rtol=residual_rtol, - residual_atol=residual_atol, - bound_tolerance=bound_tolerance, - fallback=fallback, - ) - self.config.damage_update = mode - self._last_damage_route = mode - return self - - def clear_damage_predictor(self): - """Restore the classical damage-subproblem route.""" - self.damage_controller = None - self.config.damage_update = 'classical' - self._last_damage_route = 'classical' - return self - - def _damage_step_context(self, damage_previous): - material_fields = { - name: getattr(self.material, name, None) - for name in ( - 'E', 'nu', 'Gc', 'l0', 'rho', 'eta_residual', - 'plane_stress', 'kinematics', - ) - } - time_value = ( - None if self.dt is None - else float(self._step_count * self.dt) - ) - return DamageStepContext( - step=int(self._step_count), - time=time_value, - load_factor=float(getattr(self.bcs, 'load_factor', 1.0)), - nodes=self.mesh.nodes.detach(), - elements=self.mesh.elements.detach(), - displacement=self.u.detach(), - velocity=self.v.detach(), - history_element=self.H_elem.detach(), - history_nodal=self.H_nodal.detach(), - damage_previous=damage_previous.detach(), - material=material_fields, - phase_field_model=str( - getattr(self.material, 'pf_model', 'AT2')), - energy_split=str( - getattr(self.material, 'energy_split', 'unknown')), - device=self.device, - dtype=self.dtype, - ) - def route_report(self): """Return the selected mechanics, damage, device, and physics route.""" backend_resolved = getattr(self.mechanics, 'last_backend', None) damage = { 'model': str(getattr(self.material, 'pf_model', 'AT2')), - 'update': self.config.damage_update, 'bounds': self.config.bounds_method, 'preconditioner': self.config.preconditioner, - 'last_route': self._last_damage_route, } - if self.damage_controller is not None: - damage.update(self.damage_controller.summary()) return { 'device': str(self.device), 'dtype': str(self.dtype), @@ -1587,30 +1462,19 @@ def print_route_report(self): f" fracture={report['phase_field_model']}; " f"energy_split={report['energy_split']}; " f"history_update={report['history_update']}; " - f"damage_update={damage['update']}; " f"bounds={damage['bounds']}; " f"preconditioner={damage['preconditioner']}", flush=True, ) - if self.damage_controller is not None: - print( - f" predictor={damage['predictor']}; " - f"residual_rtol={damage['residual_rtol']:.3e}; " - f"residual_atol={damage['residual_atol']:.3e}; " - f"bound_tolerance={damage['bound_tolerance']:.3e}; " - f"fallback={damage['fallback_enabled']}", - flush=True, - ) return report def get_state(self) -> dict: """Get full solver state as a dict (for checkpointing/rollback). - Hybrid solver-DL rollouts use this state for path-dependent rollback. - Keep it broader than the primary variables: rejected learned proposals - must also restore history diagnostics and installed heterogeneous - material fields so the trusted correction does not inherit stale - internal state. + The state supports path-dependent rollback. Keep it broader than the + primary variables: a restored step must also recover the history + diagnostics and installed heterogeneous material fields so that the + repeated solve does not inherit stale internal state. """ state = { 'u': self.u.clone(), @@ -2121,71 +1985,3 @@ def run(self, num_steps: int = None, output_dir: str = None, print("Done.") self.steps = int(num_steps) return history - - -class AIBridge: - """Compatibility adapter for one externally supplied damage prediction. - - New integrations should use :meth:`StaggeredSolver.set_damage_predictor`. - This adapter now applies the same admissibility and projected-residual - audit as the public learned-replacement route. - - Parameters - ---------- - solver : StaggeredSolver - residual_threshold : float - Maximum relative projected residual. Default 1e-3. - """ - - def __init__(self, solver, residual_threshold=1e-3): - self.solver = solver - self.threshold = residual_threshold - self.nn_steps = 0 - self.fallback_steps = 0 - - def step_with_nn(self, d_nn): - """Accept neural operator damage prediction, validate, maybe correct. - - Parameters - ---------- - d_nn : (N,) damage prediction from neural operator - - Returns - ------- - d_final : (N,) validated damage field - used_solver : bool — True if FEM solver was needed - """ - controller = DamageUpdateController( - lambda _context: d_nn, - mode='learned_replacement', - residual_rtol=self.threshold, - fallback=True, - ) - pf_mask = pf_vals = None - if getattr(self.solver.bcs, 'pf_dirichlet_bcs', None): - pf_mask, pf_vals = ( - self.solver.bcs.get_pf_dirichlet_mask_values()) - decision = controller.decide( - self.solver._damage_step_context(self.solver.d), - damage_solver=self.solver.damage_solver, - phase_field_mask=pf_mask, - phase_field_values=pf_vals, - ) - if decision.accepted_replacement: - self.solver.d = decision.candidate - self.solver._apply_pf_dirichlet() - self.nn_steps += 1 - return self.solver.d, False - self.solver.d = self.solver.damage_solver.solve( - self.solver.H_elem, self.solver.d) - self.fallback_steps += 1 - return self.solver.d, True - - def stats(self): - total = self.nn_steps + self.fallback_steps - if total == 0: - return "AIBridge: no steps taken" - pct = self.nn_steps / total * 100 - return (f"AIBridge: {total} steps, " - f"{self.nn_steps} NN ({pct:.0f}%), " - f"{self.fallback_steps} fallback ({100-pct:.0f}%)") diff --git a/src/phast/training/__init__.py b/src/phast/training/__init__.py deleted file mode 100644 index 0a61e512a..000000000 --- a/src/phast/training/__init__.py +++ /dev/null @@ -1,35 +0,0 @@ -"""Training utilities (standalone scaffolds, not yet wired into a trainer).""" - -from phast.training.losses import ( - adaptive_loss_weights, - dice_loss, - dice_plus_focal, - focal_loss, - log_l2_loss, - log_relative_error_loss, - rprop_step_factor, - weighted_loss_sum, -) -from phast.training.self_training import ( - SelfTrainingConfig, - filter_by_quality, - inverse_age_weights, - merge_with_labelled, - pseudo_label_step, -) - -__all__ = [ - "adaptive_loss_weights", - "dice_loss", - "dice_plus_focal", - "focal_loss", - "log_l2_loss", - "log_relative_error_loss", - "rprop_step_factor", - "weighted_loss_sum", - "SelfTrainingConfig", - "filter_by_quality", - "inverse_age_weights", - "merge_with_labelled", - "pseudo_label_step", -] diff --git a/src/phast/training/losses.py b/src/phast/training/losses.py deleted file mode 100644 index 7d3b0b873..000000000 --- a/src/phast/training/losses.py +++ /dev/null @@ -1,161 +0,0 @@ -"""Advanced loss functions for damage prediction (#51). - -Standalone module — not yet wired into a training pipeline. The training -pipeline that consumes these is separate work. - -Implements literature quick-wins from issue #51: -- log(E) loss (Manav 2024) -> log_l2_loss, log_relative_error_loss -- adaptive loss balancing (Kiyani 25) -> adaptive_loss_weights, weighted_loss_sum -- RPROP step factor (Manav 2024) -> rprop_step_factor -- Dice + Focal (Hamdi 2026) -> dice_loss, focal_loss, dice_plus_focal - -All functions are float64-compatible and autograd-friendly. -""" - -from __future__ import annotations - -import torch - - -def log_l2_loss( - pred: torch.Tensor, target: torch.Tensor, eps: float = 1e-8 -) -> torch.Tensor: - """log(L2) loss — bounds gradient magnitude across orders of magnitude. - - Returns: ``log(eps + mean((pred - target)**2))``. - """ - return torch.log(eps + ((pred - target) ** 2).mean()) - - -def log_relative_error_loss( - pred: torch.Tensor, target: torch.Tensor, eps: float = 1e-8 -) -> torch.Tensor: - """log relative error — invariant to absolute scale of target. - - Returns: ``log(eps + mean(|pred - target|) / (mean(|target|) + eps))``. - """ - num = (pred - target).abs().mean() - den = target.abs().mean() + eps - return torch.log(eps + num / den) - - -def dice_loss( - pred: torch.Tensor, target: torch.Tensor, smooth: float = 1.0 -) -> torch.Tensor: - """Dice loss for damage segmentation: ``1 - 2|P intersect T| / (|P| + |T|)``. - - Both ``pred`` and ``target`` should be in [0, 1]. Treats damage as a soft mask. - """ - p = pred.flatten() - t = target.flatten() - intersection = (p * t).sum() - return 1.0 - (2.0 * intersection + smooth) / (p.sum() + t.sum() + smooth) - - -def focal_loss( - pred: torch.Tensor, - target: torch.Tensor, - alpha: float = 0.25, - gamma: float = 2.0, - eps: float = 1e-8, -) -> torch.Tensor: - """Binary focal loss — down-weights well-classified examples. - - Both inputs assumed to be probabilities in [0, 1]. Useful for crack - localisation where most elements have ``d ~ 0`` (easy negatives). - """ - pred = pred.clamp(eps, 1.0 - eps) - is_positive = target > 0.5 - pt = torch.where(is_positive, pred, 1.0 - pred) - alpha_t = torch.where( - is_positive, - torch.as_tensor(alpha, dtype=pred.dtype, device=pred.device), - torch.as_tensor(1.0 - alpha, dtype=pred.dtype, device=pred.device), - ) - return -(alpha_t * (1.0 - pt) ** gamma * pt.log()).mean() - - -def dice_plus_focal( - pred: torch.Tensor, - target: torch.Tensor, - lambda_dice: float = 1.0, - lambda_focal: float = 1.0, - alpha: float = 0.25, - gamma: float = 2.0, -) -> torch.Tensor: - """Combined Dice + Focal — popular for class-imbalanced segmentation.""" - return lambda_dice * dice_loss(pred, target) + lambda_focal * focal_loss( - pred, target, alpha, gamma - ) - - -def rprop_step_factor( - grad: torch.Tensor, - prev_grad: torch.Tensor, - eta_plus: float = 1.2, - eta_minus: float = 0.5, - step_min: float = 1e-6, - step_max: float = 50.0, -) -> torch.Tensor: - """RPROP-style adaptive step factor. - - For each parameter, multiply its step size by ``eta_plus`` when ``grad`` - sign matches ``prev_grad``, ``eta_minus`` otherwise. Returns the - multiplicative step factor, clamped to ``[step_min, step_max]``. - """ - sign_change = grad * prev_grad - factor = torch.where( - sign_change > 0, - torch.full_like(grad, eta_plus), - torch.where( - sign_change < 0, - torch.full_like(grad, eta_minus), - torch.ones_like(grad), - ), - ) - return factor.clamp(min=step_min, max=step_max) - - -def adaptive_loss_weights( - losses: dict[str, torch.Tensor], eps: float = 1e-8 -) -> dict[str, torch.Tensor]: - """Inverse-magnitude weighting (Kiyani 2025). - - Returns detached tensor weights ``w_i = (1 / (|L_i| + eps)) / sum_j(...)`` - on the same dtype/device as the first loss. The weights do not propagate - gradients, but they remain tensors so callers can use them directly in a - torch training loop without losing dtype/device information. - """ - if not losses: - raise ValueError("losses must contain at least one term") - first = next(iter(losses.values())) - dtype = first.dtype - device = first.device - inv_terms: dict[str, torch.Tensor] = {} - with torch.no_grad(): - for name, loss in losses.items(): - mag = loss.detach().abs().to(dtype=dtype, device=device) - inv_terms[name] = 1.0 / (mag + eps) - denom = torch.stack([v.reshape(()) for v in inv_terms.values()]).sum() - weights = {name: inv / denom for name, inv in inv_terms.items()} - return weights - - -def weighted_loss_sum( - losses: dict[str, torch.Tensor], - weights: dict[str, torch.Tensor] | None = None, - eps: float = 1e-8, -) -> torch.Tensor: - """Return ``sum_i w_i L_i`` using adaptive weights by default.""" - if not losses: - raise ValueError("losses must contain at least one term") - if weights is None: - weights = adaptive_loss_weights(losses, eps=eps) - total = None - for name, loss in losses.items(): - if name not in weights: - raise KeyError(f"missing weight for loss term {name!r}") - term = weights[name].to(dtype=loss.dtype, device=loss.device) * loss - total = term if total is None else total + term - assert total is not None - return total diff --git a/src/phast/training/self_training.py b/src/phast/training/self_training.py deleted file mode 100644 index cc811c373..000000000 --- a/src/phast/training/self_training.py +++ /dev/null @@ -1,177 +0,0 @@ -"""Self-training augmentation for damage rollout prediction (#56). - -Pipeline: after the teacher-forced phase, run the model autoregressively, -filter rollout snapshots by a quality criterion (e.g. PDE residual below -threshold or confidence above threshold), and append them to the training -set with an inverse-age decay weight ``w = 1 / (1 + age)``. - -References ----------- -- Yarowsky (1995) — bootstrapping for word-sense disambiguation -- Lee (2013) — pseudo-label deep learning -- Sohn et al. (2020) — FixMatch (consistency + confidence filtering) -- Schwarzer et al. (2019) — temporal downsampling + self-training (rollout - predictions added back into the training set; the source for #56) - -Standalone scaffold — not wired into a training loop. All ops are -float64-compatible and autograd-friendly where applicable; the filtering -step is wrapped in ``torch.no_grad()`` because pseudo-labels enter the -graph as constants (the gradient flow during retraining comes from the -standard supervised loss, not from the labelling step). -""" - -from __future__ import annotations - -from dataclasses import dataclass -from typing import Callable, Optional, Tuple - -import torch - - -@dataclass -class SelfTrainingConfig: - """Knobs for the self-training augmentation pass.""" - - quality_threshold: float = 0.9 - """For confidence mode: minimum predicted confidence in [0, 1]. - For residual mode: maximum allowed PDE residual (inputs above are dropped). - """ - use_residual_filter: bool = False - """If True, filter by a user-supplied PDE-residual callable (Schwarzer 2019); - else filter by output confidence (Lee 2013). - """ - use_entropy_confidence: bool = True - """When in confidence mode and outputs are class-prob: if True use - ``1 - entropy / log(C)``, else use max-prob. Ignored for scalar outputs. - """ - max_pseudo_per_pass: int = 10000 - """Cap on number of pseudo-labelled samples appended per pass.""" - age_decay: bool = True - """If True, weight pseudo samples by ``1 / (1 + age)`` (Schwarzer 2019).""" - - -def _confidence_from_probs( - probs: torch.Tensor, use_entropy: bool -) -> torch.Tensor: - """(n, C) probabilities -> (n,) confidence scalar in [0, 1].""" - if use_entropy: - eps = torch.tensor(1e-12, dtype=probs.dtype, device=probs.device) - clamped = probs.clamp(min=eps) - entropy = -(clamped * clamped.log()).sum(dim=-1) - n_classes = probs.shape[-1] - log_c = torch.log( - torch.tensor(float(n_classes), dtype=probs.dtype, device=probs.device) - ) - return 1.0 - entropy / log_c - return probs.max(dim=-1).values - - -def filter_by_quality( - scores: torch.Tensor, - config: SelfTrainingConfig, -) -> torch.Tensor: - """Boolean mask (n,) selecting samples that pass the quality filter. - - ``scores`` is *confidence* (higher better, threshold is lower bound) when - ``use_residual_filter=False``, else *residual* (lower better, threshold - is upper bound). - """ - if config.use_residual_filter: - return scores <= config.quality_threshold - return scores >= config.quality_threshold - - -def inverse_age_weights( - ages: torch.Tensor, -) -> torch.Tensor: - """Inverse-age decay weights ``w = 1 / (1 + age)`` (Schwarzer 2019).""" - return 1.0 / (1.0 + ages.to(dtype=torch.float64)) - - -def pseudo_label_step( - model_predict: Callable[[torch.Tensor], torch.Tensor], - unlabelled_inputs: torch.Tensor, - config: SelfTrainingConfig, - residual_fn: Optional[Callable[[torch.Tensor, torch.Tensor], torch.Tensor]] = None, -) -> Tuple[torch.Tensor, torch.Tensor, torch.Tensor]: - """Run one pseudo-labelling pass. - - Parameters - ---------- - model_predict - Callable returning model output for a batch of inputs. May return - scalar regression output of shape (n, ...) or class probabilities - of shape (n, C). - unlabelled_inputs - (n, ...) batch of unlabelled features (e.g. autoregressive rollout - snapshots that were not seen during teacher-forced training). - config - Filter knobs. - residual_fn - Required when ``config.use_residual_filter=True``. Maps - ``(inputs, predictions) -> (n,)`` per-sample residual scalars. - - Returns - ------- - (filtered_inputs, pseudo_labels, scores) - with ``n_kept <= n``. ``scores`` are the confidence/residual values - used for filtering; useful for downstream weighting. - """ - with torch.no_grad(): - preds = model_predict(unlabelled_inputs) - if config.use_residual_filter: - if residual_fn is None: - raise ValueError("residual_fn required when use_residual_filter=True") - scores = residual_fn(unlabelled_inputs, preds) - else: - if preds.dim() >= 2 and preds.shape[-1] > 1: - scores = _confidence_from_probs(preds, config.use_entropy_confidence) - else: - # Scalar regression: treat distance from {0, 1} as a - # binary-confidence proxy (damage ∈ [0, 1] in PF problems). - p = preds.reshape(preds.shape[0], -1).mean(dim=-1) - scores = (p - 0.5).abs() * 2.0 - mask = filter_by_quality(scores, config) - if int(mask.sum().item()) > config.max_pseudo_per_pass: - kept_idx = torch.nonzero(mask, as_tuple=True)[0] - if config.use_residual_filter: - topk = torch.topk( - -scores[kept_idx], k=config.max_pseudo_per_pass - ).indices - else: - topk = torch.topk( - scores[kept_idx], k=config.max_pseudo_per_pass - ).indices - new_mask = torch.zeros_like(mask) - new_mask[kept_idx[topk]] = True - mask = new_mask - if preds.dim() >= 2 and preds.shape[-1] > 1 and not config.use_residual_filter: - pseudo_labels = preds.argmax(dim=-1) - else: - pseudo_labels = preds - return unlabelled_inputs[mask], pseudo_labels[mask], scores[mask] - - -def merge_with_labelled( - labelled_x: torch.Tensor, - labelled_y: torch.Tensor, - pseudo_x: torch.Tensor, - pseudo_y: torch.Tensor, - pseudo_ages: Optional[torch.Tensor] = None, -) -> Tuple[torch.Tensor, torch.Tensor, torch.Tensor]: - """Concatenate labelled + pseudo into a single combined dataset. - - Pseudo samples are appended last so their per-sample weights can be - addressed via slicing. Real samples receive weight 1.0; pseudo samples - receive ``1 / (1 + age)`` if ``pseudo_ages`` is provided, else 1.0. - """ - x = torch.cat([labelled_x, pseudo_x], dim=0) - y = torch.cat([labelled_y, pseudo_y], dim=0) - real_w = torch.ones(labelled_x.shape[0], dtype=torch.float64) - pseudo_w = ( - inverse_age_weights(pseudo_ages) - if pseudo_ages is not None - else torch.ones(pseudo_x.shape[0], dtype=torch.float64) - ) - weights = torch.cat([real_w, pseudo_w], dim=0) - return x, y, weights diff --git a/src/phast/utils/doctor.py b/src/phast/utils/doctor.py index d71e91610..016282495 100644 --- a/src/phast/utils/doctor.py +++ b/src/phast/utils/doctor.py @@ -106,8 +106,6 @@ def build_report() -> str: " spectral/Amor QS on CPU/HPC: install PETSc/MUMPS so backend=auto can select mumps", " cohesive contact: sparse quasi-static backend, backend=auto, normal-contact penalty only when configured", " J2 plasticity: sparse quasi-static backend, backend=auto, guarded supported material combinations", - " dataset/deep learning: Zarr trajectory stores; MP4/raster visualisation when animations are requested", - " learned damage: classical is default; learned_proposal keeps exact correction; learned_replacement is audited experimental inference", "", "Install guidance:", " Start with: pip install -e .", diff --git a/src/phast/utils/io_utils.py b/src/phast/utils/io_utils.py index a9d4d88dd..9aea49c9f 100644 --- a/src/phast/utils/io_utils.py +++ b/src/phast/utils/io_utils.py @@ -2,8 +2,8 @@ Output utilities: HDF5 trajectories, VTU (ParaView), CSV (history), and explicitly selected Zarr trajectory stores. -Includes helpers for GNO (Graph Neural Operator) training data export: -- ``compute_edge_index`` converts triangle connectivity to PyG edge_index. +Includes mesh-graph and run-bookkeeping helpers: +- ``compute_edge_index`` converts triangle connectivity to a COO edge list. - ``write_profiler_csv`` exports solver profiling data. - ``generate_run_tag`` creates informative output directory names. - ``save_run_metadata`` writes a JSON with full run context. @@ -127,7 +127,7 @@ def plot_energy_history(rows, path: str, xlabel: str = "Step") -> None: def compute_edge_index(elements): - """Compute PyG-style edge_index from triangle connectivity. + """Compute a ``(2, n_edges)`` COO edge index from triangle connectivity. For each triangle (i, j, k), creates bidirectional edges: (i,j), (j,i), (j,k), (k,j), (i,k), (k,i). @@ -777,7 +777,7 @@ def _get_step_array(name): def init_h5(h5_path: str, mesh, material): """Create and initialize an H5 file with mesh and metadata. - Stores node coordinates, element connectivity, PyG-compatible + Stores node coordinates, element connectivity, the COO-format mesh ``edge_index``, and mesh size attributes (``n_nodes``, ``n_elements``). Returns @@ -796,7 +796,7 @@ def init_h5(h5_path: str, mesh, material): mesh_grp.create_dataset('element_connectivity', data=mesh.elements.cpu().numpy()) - # PyG-compatible edge_index for GNO training + # Undirected mesh edges in COO (2, n_edges) layout edge_index = compute_edge_index(mesh.elements) mesh_grp.create_dataset('edge_index', data=edge_index) diff --git a/src/phast/workflow/execution.py b/src/phast/workflow/execution.py index d0e50ef42..a96d1799b 100644 --- a/src/phast/workflow/execution.py +++ b/src/phast/workflow/execution.py @@ -558,10 +558,6 @@ def _fracture_bridge_issues(spec: ProblemSpec) -> tuple[str, ...]: solver = spec.solver.parameters defaults = LegacySolverSettings(solver_type=spec.solver.kind) - if solver.get("damage_update", "classical") != "classical" or any( - solver.get(key) for key in ("damage_predictor", "damage_checkpoint", "damage_predictor_options") - ): - issues.append("learned damage execution is not supported by this schema-v2 adapter") if solver.get("adaptive_dt", False): issues.append("adaptive_dt is unsupported by the fixed-time CLI loading/output schedule") try: diff --git a/tests/test_damage_step_context_fem_helpers.py b/tests/test_damage_step_context_fem_helpers.py deleted file mode 100644 index 62dcee382..000000000 --- a/tests/test_damage_step_context_fem_helpers.py +++ /dev/null @@ -1,122 +0,0 @@ -"""Tests for the architecture-neutral finite-element helpers on the context. - -These helpers are shared by every learned-damage predictor, whatever its -architecture, so they are tested here rather than beside any one adapter. -""" -from __future__ import annotations - -import pytest -import torch - -from phast.learned_damage import DamageStepContext - -DTYPE = torch.float64 - -# Unit square split into four triangles around a central node. -NODES = torch.tensor( - [[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0], [0.5, 0.5]], dtype=DTYPE) -ELEMENTS = torch.tensor( - [[0, 1, 4], [1, 2, 4], [2, 3, 4], [3, 0, 4]], dtype=torch.long) - - -def make_context(nodes=NODES, elements=ELEMENTS, *, history=1.0): - n_nodes = nodes.shape[0] - return DamageStepContext( - step=1, - time=0.0, - load_factor=1.0, - nodes=nodes, - elements=elements, - displacement=torch.zeros((n_nodes, 2), dtype=DTYPE), - velocity=torch.zeros((n_nodes, 2), dtype=DTYPE), - history_element=torch.full((elements.shape[0],), history, dtype=DTYPE), - history_nodal=torch.zeros((n_nodes,), dtype=DTYPE), - damage_previous=torch.zeros((n_nodes,), dtype=DTYPE), - material={"Gc": 2.7, "l0": 0.4, "pf_model": "AT2"}, - phase_field_model="AT2", - energy_split="spectral", - device=torch.device("cpu"), - dtype=DTYPE, - ) - - -def test_element_areas_sum_to_the_domain_area(): - assert float(make_context().element_areas().sum()) == pytest.approx(1.0) - - -def test_element_areas_are_orientation_independent(): - flipped = ELEMENTS[:, [0, 2, 1]] - assert torch.allclose( - make_context().element_areas(), - make_context(elements=flipped).element_areas()) - - -def test_assembly_conserves_the_integral(): - context = make_context() - values = torch.tensor([1.0, 2.0, 3.0, 4.0], dtype=DTYPE) - assembled = context.assemble_element_to_nodes(values) - expected = float((values * context.element_areas()).sum()) - assert float(assembled.sum()) == pytest.approx(expected) - - -def test_assembly_of_a_unit_field_returns_the_lumped_mass(): - context = make_context() - lumped = context.assemble_element_to_nodes( - torch.ones(ELEMENTS.shape[0], dtype=DTYPE)) - # The central node belongs to every triangle, so it carries a third of - # the whole domain area. - assert float(lumped[4]) == pytest.approx(1.0 / 3.0) - assert float(lumped.sum()) == pytest.approx(1.0) - - -def test_assembly_rejects_a_wrong_length(): - with pytest.raises(ValueError, match="one entry per element"): - make_context().assemble_element_to_nodes(torch.zeros(2, dtype=DTYPE)) - - -def test_boundary_mask_marks_the_perimeter_and_not_the_interior(): - assert make_context().boundary_node_mask().tolist() == [ - True, True, True, True, False] - - -def test_boundary_mask_detects_an_interior_void(): - # An annulus-like patch: a square ring of eight triangles around a hole. - outer = torch.tensor( - [[0.0, 0.0], [1.0, 0.0], [2.0, 0.0], [2.0, 1.0], - [2.0, 2.0], [1.0, 2.0], [0.0, 2.0], [0.0, 1.0]], dtype=DTYPE) - inner = torch.tensor( - [[0.75, 0.75], [1.25, 0.75], [1.25, 1.25], [0.75, 1.25]], dtype=DTYPE) - nodes = torch.cat((outer, inner), dim=0) - ring = torch.tensor( - [[0, 1, 8], [1, 2, 9], [1, 9, 8], [2, 3, 9], - [3, 4, 10], [3, 10, 9], [4, 5, 10], [5, 6, 11], - [5, 11, 10], [6, 7, 11], [7, 0, 8], [7, 8, 11]], dtype=torch.long) - mask = make_context(nodes, ring).boundary_node_mask() - # Every node lies on either the outer perimeter or the hole. - assert bool(mask.all()) - - -def test_edge_lengths_match_the_geometry(): - context = make_context() - edge_index = torch.tensor([[0, 0], [1, 4]], dtype=torch.long) - lengths = context.edge_lengths(edge_index) - assert float(lengths[0]) == pytest.approx(1.0) - assert float(lengths[1]) == pytest.approx(0.5 * 2.0 ** 0.5) - - -def test_graph_edge_index_is_symmetric_and_duplicate_free(): - edge_index = make_context().graph_edge_index() - pairs = {tuple(pair) for pair in edge_index.t().tolist()} - assert len(pairs) == edge_index.shape[1] - assert all((b, a) in pairs for a, b in pairs) - - -def test_canonical_node_features_have_six_columns(): - features = make_context().canonical_node_features() - assert features.shape == (NODES.shape[0], 6) - - -def test_element_areas_reject_a_non_triangular_mesh(): - quads = torch.tensor([[0, 1, 2, 3]], dtype=torch.long) - with pytest.raises(NotImplementedError, match="three-node triangles"): - make_context(elements=quads).element_areas() diff --git a/tests/test_learned_damage.py b/tests/test_learned_damage.py deleted file mode 100644 index 3611a0e37..000000000 --- a/tests/test_learned_damage.py +++ /dev/null @@ -1,152 +0,0 @@ -from __future__ import annotations - -import pytest -import torch - -from phast.learned_damage import ( - DamagePrediction, - DamagePredictionRejected, - DamageStepContext, - DamageUpdateController, -) - - -class ConstantPredictor: - name = "constant" - - def __init__(self, value): - self.value = value - - def predict(self, context): - return DamagePrediction( - damage=torch.full_like(context.damage_previous, self.value) - ) - - -class IncrementPredictor: - name = "increment" - - def predict(self, context): - return DamagePrediction( - damage=torch.full_like(context.damage_previous, 0.2), - representation="damage_increment", - ) - - -class ResidualSolver: - def __init__(self, target): - self.target = target - - def compute_residual(self, history, damage): - del history - return damage - self.target - - -def make_context(previous=0.2): - dtype = torch.float64 - previous_damage = torch.full((3,), previous, dtype=dtype) - return DamageStepContext( - step=2, - time=0.1, - load_factor=0.5, - nodes=torch.tensor([[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]], dtype=dtype), - elements=torch.tensor([[0, 1, 2]], dtype=torch.long), - displacement=torch.zeros((3, 2), dtype=dtype), - velocity=torch.zeros((3, 2), dtype=dtype), - history_element=torch.zeros((1,), dtype=dtype), - history_nodal=torch.zeros((3,), dtype=dtype), - damage_previous=previous_damage, - material={"pf_model": "AT2", "energy_split": "spectral"}, - phase_field_model="AT2", - energy_split="spectral", - device=torch.device("cpu"), - dtype=dtype, - ) - - -def make_controller(mode, value, target, fallback=True): - controller = DamageUpdateController( - predictor=ConstantPredictor(value), - mode=mode, - fallback=fallback, - residual_rtol=1.0e-6, - residual_atol=1.0e-12, - bound_tolerance=1.0e-8, - ) - solver = ResidualSolver(torch.full((3,), target, dtype=torch.float64)) - return controller, solver - - -def test_proposal_is_projected_and_returned_as_initial_guess(): - context = make_context(previous=0.2) - controller, damage_solver = make_controller( - mode="learned_proposal", value=0.1, target=0.4 - ) - decision = controller.decide(context, damage_solver=damage_solver) - - assert decision.route == "learned_proposal_exact_correction" - assert not decision.accepted_replacement - assert torch.equal(decision.candidate, context.damage_previous) - - -def test_replacement_accepts_audited_prediction(): - context = make_context(previous=0.2) - controller, damage_solver = make_controller( - mode="learned_replacement", value=0.4, target=0.4 - ) - decision = controller.decide(context, damage_solver=damage_solver) - - assert decision.route == "learned_replacement" - assert decision.accepted_replacement - assert torch.allclose( - decision.candidate, - torch.full((3,), 0.4, dtype=torch.float64), - ) - - -def test_replacement_accepts_damage_increment_representation(): - context = make_context(previous=0.2) - controller = DamageUpdateController( - predictor=IncrementPredictor(), - mode="learned_replacement", - fallback=True, - residual_rtol=1.0e-6, - residual_atol=1.0e-12, - bound_tolerance=1.0e-8, - ) - damage_solver = ResidualSolver( - torch.full((3,), 0.4, dtype=torch.float64) - ) - - decision = controller.decide(context, damage_solver=damage_solver) - assert decision.accepted_replacement - assert torch.allclose( - decision.candidate, - torch.full((3,), 0.4, dtype=torch.float64), - ) - - -def test_replacement_rejects_invalid_prediction_and_requests_fallback(): - context = make_context(previous=0.2) - controller, damage_solver = make_controller( - mode="learned_replacement", value=float("nan"), target=0.4 - ) - decision = controller.decide(context, damage_solver=damage_solver) - - assert decision.route == "classical_fallback" - assert not decision.accepted_replacement - assert decision.candidate is None - assert "NaN" in decision.reason - - -def test_replacement_can_raise_instead_of_falling_back(): - context = make_context(previous=0.2) - controller, damage_solver = make_controller( - mode="learned_replacement", - value=float("nan"), - target=0.4, - fallback=False, - ) - - with pytest.raises(DamagePredictionRejected): - controller.decide(context, damage_solver=damage_solver) diff --git a/tests/test_mesh_graph_net_adapter.py b/tests/test_mesh_graph_net_adapter.py deleted file mode 100644 index 170cb3a35..000000000 --- a/tests/test_mesh_graph_net_adapter.py +++ /dev/null @@ -1,251 +0,0 @@ -"""Tests for the mesh-graph-network learned-damage adapter. - -The adapter is one architecture among several. What is checked here is the -wrapper: its feature contract, its checkpoint handling and its conformance to -the predictor protocol. The finite-element helpers it builds on are covered by -``test_damage_step_context_fem_helpers``. -""" -from __future__ import annotations - -import pytest -import torch - -from phast.learned_damage import DamageStepContext - -from examples.learned_damage.architectures.mesh_graph_net import ( - MeshGraphDamageNet, - MeshGraphNetPredictor, - bounding_box_indicator, - create_predictor, - export_torchscript, - load_network, -) - -DTYPE = torch.float64 - -NODES = torch.tensor( - [[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0], [0.5, 0.5]], dtype=DTYPE) -ELEMENTS = torch.tensor( - [[0, 1, 4], [1, 2, 4], [2, 3, 4], [3, 0, 4]], dtype=torch.long) - - -def make_context(*, history=1.0, Gc=2.7, l0=0.4): - n_nodes = NODES.shape[0] - return DamageStepContext( - step=1, - time=0.0, - load_factor=1.0, - nodes=NODES, - elements=ELEMENTS, - displacement=torch.zeros((n_nodes, 2), dtype=DTYPE), - velocity=torch.zeros((n_nodes, 2), dtype=DTYPE), - history_element=torch.full((ELEMENTS.shape[0],), history, dtype=DTYPE), - history_nodal=torch.zeros((n_nodes,), dtype=DTYPE), - damage_previous=torch.zeros((n_nodes,), dtype=DTYPE), - material={"Gc": Gc, "l0": l0, "pf_model": "AT2"}, - phase_field_model="AT2", - energy_split="spectral", - device=torch.device("cpu"), - dtype=DTYPE, - ) - - -def make_net(): - torch.manual_seed(0) - return MeshGraphDamageNet( - latent_dim=4, hidden_layers=1, message_passing_steps=2) - - -def random_graph_inputs(context): - edge_index = context.graph_edge_index() - return ( - torch.rand(NODES.shape[0], 2), - edge_index, - torch.rand(edge_index.shape[1], 1), - ) - - -def save_state_dict(path, net): - torch.save( - { - "model_state_dict": net.state_dict(), - "configuration": { - "latent_dim": 4, "hidden_layers": 1, "message_passing_steps": 2}, - "completed_epochs": 3, - }, - path, - ) - return path - - -# ---------------------------------------------------------------- network - -def test_network_returns_one_value_per_node(): - net = make_net() - with torch.no_grad(): - output = net(*random_graph_inputs(make_context())) - assert output.shape == (NODES.shape[0],) - - -def test_network_is_scriptable_and_agrees_with_the_eager_model(): - net = make_net().eval() - inputs = random_graph_inputs(make_context()) - scripted = torch.jit.script(net) - with torch.no_grad(): - assert torch.equal(net(*inputs), scripted(*inputs)) - - -# ---------------------------------------------------------------- features - -def test_bounding_box_indicator_ignores_the_interior(): - assert bounding_box_indicator(NODES).tolist() == [1.0, 1.0, 1.0, 1.0, 0.0] - - -def test_node_features_use_the_assembled_normalized_history(): - predictor = MeshGraphNetPredictor(make_net()) - context = make_context(history=1.0, Gc=2.7, l0=0.4) - - features = predictor.node_features(context) - expected = context.assemble_element_to_nodes( - context.history_element / (2.7 * 0.4)) - - assert features.shape == (NODES.shape[0], 2) - assert torch.allclose(features[:, 0], expected) - assert torch.allclose(features[:, 1], bounding_box_indicator(NODES)) - - -def test_topological_boundary_indicator_can_be_selected(): - context = make_context() - predictor = MeshGraphNetPredictor( - make_net(), boundary_indicator="topological") - features = predictor.node_features(context) - assert torch.allclose( - features[:, 1], context.boundary_node_mask().to(DTYPE)) - - -def test_an_unknown_boundary_indicator_is_rejected(): - with pytest.raises(ValueError, match="bounding_box.*topological"): - MeshGraphNetPredictor(make_net(), boundary_indicator="guess") - - -def test_overrides_take_precedence_over_the_running_material(): - context = make_context(Gc=2.7, l0=0.4) - default = MeshGraphNetPredictor(make_net()).node_features(context) - overridden = MeshGraphNetPredictor( - make_net(), Gc=1.0, length_scale=1.0).node_features(context) - assert not torch.allclose(default[:, 0], overridden[:, 0]) - - -def test_a_material_without_the_normalization_constants_is_rejected(): - stripped = DamageStepContext( - **{**make_context().__dict__, "material": {"pf_model": "AT2"}}) - with pytest.raises(ValueError, match="exposes no 'Gc'"): - MeshGraphNetPredictor(make_net()).node_features(stripped) - - -def test_graph_is_cached_between_calls_on_one_mesh(): - predictor = MeshGraphNetPredictor(make_net()) - context = make_context() - assert predictor._graph(context)[0] is predictor._graph(context)[0] - - -def test_graph_cache_invalidates_for_a_different_same_sized_mesh(): - predictor = MeshGraphNetPredictor(make_net()) - context = make_context() - first_lengths = predictor._graph(context)[1] - scaled = DamageStepContext( - **{**context.__dict__, "nodes": context.nodes * 2.0}) - second_lengths = predictor._graph(scaled)[1] - assert not torch.equal(first_lengths, second_lengths) - - -# ---------------------------------------------------------------- protocol - -def test_predict_returns_one_nodal_value_in_the_context_dtype(): - context = make_context() - prediction = MeshGraphNetPredictor(make_net()).predict(context) - assert prediction.damage.shape == context.damage_previous.shape - assert prediction.damage.dtype == context.dtype - assert prediction.diagnostics["adapter"] == "mesh-graph-net" - assert prediction.diagnostics["boundary_indicator"] == "bounding_box" - - -def test_increment_representation_is_passed_through(): - predictor = MeshGraphNetPredictor( - make_net(), representation="damage_increment") - assert predictor.predict(make_context()).representation == "damage_increment" - - -# ---------------------------------------------------------------- checkpoints - -def test_state_dict_checkpoint_round_trips(tmp_path): - net = make_net().eval() - loaded, metadata = load_network(save_state_dict(tmp_path / "ck.pt", net)) - - assert metadata["format"] == "state_dict" - assert metadata["completed_epochs"] == 3 - - inputs = random_graph_inputs(make_context()) - with torch.no_grad(): - assert torch.equal(net(*inputs), loaded(*inputs)) - - -def test_a_checkpoint_without_architecture_entries_is_rejected(tmp_path): - path = tmp_path / "incomplete.pt" - torch.save({"model_state_dict": make_net().state_dict()}, path) - with pytest.raises(ValueError, match="missing architecture entries"): - load_network(path) - - -def test_an_unrecognized_file_is_rejected(tmp_path): - path = tmp_path / "unrelated.pt" - torch.save({"something": 1}, path) - with pytest.raises(ValueError, match="neither a TorchScript archive"): - load_network(path) - - -def test_torchscript_export_loads_and_predicts(tmp_path): - net = make_net().eval() - archive = export_torchscript( - save_state_dict(tmp_path / "ck.pt", net), - tmp_path / "scripted.pt", - n_nodes=16, - ) - - loaded, metadata = load_network(archive) - assert metadata["format"] == "torchscript" - - inputs = random_graph_inputs(make_context()) - with torch.no_grad(): - assert torch.equal(net(*inputs), loaded(*inputs)) - - predictor = create_predictor(checkpoint=archive, device="cpu", options={}) - assert predictor.predict(make_context()).damage.shape == (NODES.shape[0],) - - -def test_exporting_an_archive_again_is_rejected(tmp_path): - archive = export_torchscript( - save_state_dict(tmp_path / "ck.pt", make_net().eval()), - tmp_path / "scripted.pt", - n_nodes=16, - ) - with pytest.raises(ValueError, match="already a TorchScript archive"): - export_torchscript(archive, tmp_path / "again.pt") - - -def test_factory_passes_options_through(tmp_path): - archive = export_torchscript( - save_state_dict(tmp_path / "ck.pt", make_net().eval()), - tmp_path / "scripted.pt", - n_nodes=16, - ) - predictor = create_predictor( - checkpoint=archive, - device="cpu", - options={ - "representation": "damage_increment", - "boundary_indicator": "topological", - }, - ) - assert predictor.representation == "damage_increment" - assert predictor.boundary_indicator == "topological" diff --git a/tests/test_standard_workflow.py b/tests/test_standard_workflow.py index 3f35711be..1fff9816d 100644 --- a/tests/test_standard_workflow.py +++ b/tests/test_standard_workflow.py @@ -161,7 +161,7 @@ def test_automatic_cfl_loading_and_hdf5_default(tmp_path: Path) -> None: @pytest.mark.parametrize("case", [ "multiple_steps", "step_mismatch", "multimaterial_dynamic", "partial_material", "geometric_region", "initial_velocity", "initial_damage_range", "unknown_solver", - "learned", "implicit_integrator", "adaptive_dt", "unused_backend", "unknown_control", + "implicit_integrator", "adaptive_dt", "unused_backend", "unknown_control", "missing_dt", "both_time_limits", "fixed_dt_and_cfl", "unknown_geometry", "3d", "higher_order", "material_typo", "negative_density", "nonfinite_material", "unknown_output", "vtu", "field_without_trajectory", "field_cadence", @@ -185,9 +185,9 @@ def test_unsupported_choices_fail_before_execution(case: str, tmp_path: Path, mo raw["initial_conditions"][0]["field"] = "velocity" elif case == "initial_damage_range": raw["initial_conditions"][0]["value"] = 2.0 - elif case in {"unknown_solver", "learned", "implicit_integrator", "adaptive_dt", "unused_backend"}: + elif case in {"unknown_solver", "implicit_integrator", "adaptive_dt", "unused_backend"}: key, value = { - "unknown_solver": ("damage_toll", 1.0e-6), "learned": ("damage_update", "learned_proposal"), + "unknown_solver": ("damage_toll", 1.0e-6), "implicit_integrator": ("time_integrator", "generalized_alpha"), "adaptive_dt": ("adaptive_dt", True), "unused_backend": ("backend", "scipy"), }[case]