Upload folder using huggingface_hub

.gitattributes

@@ -33,3 +33,10 @@ saved_model/**/* filter=lfs diff=lfs merge=lfs -text
 *.zip filter=lfs diff=lfs merge=lfs -text
 *.zst filter=lfs diff=lfs merge=lfs -text
 *tfevents* filter=lfs diff=lfs merge=lfs -text
+results/best_model/figures/fields_t1.png filter=lfs diff=lfs merge=lfs -text
+results/best_model/figures/fields_t10.png filter=lfs diff=lfs merge=lfs -text
+results/best_model/figures/fields_t5.png filter=lfs diff=lfs merge=lfs -text
+results/best_model/figures/streamlines_t1.png filter=lfs diff=lfs merge=lfs -text
+results/best_model/figures/streamlines_t10.png filter=lfs diff=lfs merge=lfs -text
+results/best_model/figures/streamlines_t5.png filter=lfs diff=lfs merge=lfs -text
+results/best_model/figures/training_curves.png filter=lfs diff=lfs merge=lfs -text

README.md

@@ -0,0 +1,86 @@
+---
+license: other
+language:
+- en
+tags:
+- OneScience
+- cavity-flow
+- physics-informed-neural-network
+- quantum-machine-learning
+---
+
+# QCPINN-CavityFlow
+
+## Model Introduction
+
+This package reproduces the unsteady lid-driven cavity experiment in [QCPINN: Quantum-Classical Physics-Informed Neural Networks for Solving PDEs](https://arxiv.org/abs/2503.16678), Table 3. The Angle-Cascade quantum-classical PINN maps `(t, x, y)` to `(u, v, p)`. Sobol-sampled PDE and boundary points drive training; the author's reference field is used for validation and evaluation, not gradient updates.
+
+Ten independent runs of 20,000 epochs were completed on CPU quantum simulation. The included `weight/best_model.pt` is the validation-selected checkpoint from seed 6, epoch 1, not a claim of paper-level accuracy. Only the reference-data path in its checkpoint metadata was changed for portability; all nine parameter tensors are unchanged. The published checkpoint SHA256 is `498bcc0dbfc053ca4824b1b6cd2914e85536d0aed1e07ad50d5ef3b4a0b11b87`; the original training checkpoint SHA256 was `f24e5f3357b045ed2a49fdb8a377cd3668a5cbb8100ff3c13906e3bfce9aff11`.
+
+## Intended Use
+
+| Task | Description |
+| --- | --- |
+| Cavity-flow field inference | Predict velocity and pressure at spatiotemporal coordinates |
+| PINN research | Examine a quantum-classical Navier-Stokes training protocol |
+
+## Usage
+
+### Download the Model Package
+
+```bash
+hf download OneScience-Group/QCPINN-CavityFlow --local-dir ./QCPINN-CavityFlow
+cd QCPINN-CavityFlow
+```
+
+### Install Dependencies
+
+The formal run used Python 3.10, CPU PyTorch 2.6.0 and PennyLane 0.40.0. The `+cpu` PyTorch wheel requires the official CPU index:
+
+```bash
+conda create -n qcpinn-cavity python=3.10 -y
+conda activate qcpinn-cavity
+pip install -r requirements.txt --extra-index-url https://download.pytorch.org/whl/cpu
+pip install 'modelscope[datasets]'
+```
+
+### Reference Data
+
+Hugging Face hosts only the model package. Obtain the author's `cavity.mat` from [OneScience/QCPINN-Cavity on ModelScope](https://www.modelscope.cn/datasets/OneScience/QCPINN-Cavity):
+
+```bash
+modelscope download --dataset OneScience/QCPINN-Cavity --local_dir ./data
+```
+
+The expected file is `data/cavity.mat`, SHA256 `4010d55dc23aa0f8ba79b24baf75e3903044bc4333d772168121cd30bc95c3ba`. It contains `cavity_internal`, `cavity_gamma0`, and `cavity_gamma1` in MATLAB v7.3/HDF5 format. No external labeled training dataset is needed.
+
+### Train
+
+From a fresh package copy on a SLURM compute cluster, with the Python environment activated:
+
+```bash
+sbatch scripts/run_all.slurm
+```
+
+The script runs ten seeds, 20,000 epochs per seed, then reloads checkpoints for inference and figures. Adjust the cluster's time/partition settings as needed. Preserve the bundled weight and results before retraining.
+
+### Inference and Results
+
+```bash
+python scripts/inference.py --checkpoint weight/best_model.pt
+python scripts/result.py
+```
+
+The selected checkpoint's internal-field relative L2 errors for `(u, v, p)` were `(123.70%, 121.43%, 164.19%)`. Across ten validation-best checkpoints, the means were `(200.42%, 274.73%, 148.86%)`, versus `(25.73%, 38.00%, 42.65%)` in paper Table 3. The full training was completed, but the paper accuracy was not reproduced. Historical `results/best_model/metrics.json` records the original checkpoint hash; rerunning inference records the published checkpoint hash.
+
+## OneScience Repositories
+
+| Platform | Main repository | Skills repository |
+| --- | --- | --- |
+| Gitee | [OneScience](https://gitee.com/onescience-ai/onescience) | [OneSkills](https://gitee.com/onescience-ai/oneskills) |
+| GitHub | [OneScience](https://github.com/onescience-ai/OneScience) | [OneSkills](https://github.com/onescience-ai/oneskills) |
+
+## Citation and License
+
+- Paper: [QCPINN](https://arxiv.org/abs/2503.16678).
+- Reference field: [authors' QCPINN repository](https://github.com/afrah/QCPINN/tree/main/data). The repository has an MIT software license; the data is distributed separately with its own provenance information.

config.json

@@ -0,0 +1,12 @@
+{
+  "model_type": "qcpinn_cavity_angle_cascade",
+  "architectures": ["QCPINN"],
+  "framework": "pytorch",
+  "experiment": "table3_unsteady_cavity",
+  "input_dim": 3,
+  "output_dim": 3,
+  "hidden_dim": 50,
+  "n_qubits": 5,
+  "checkpoint": "weight/best_model.pt",
+  "reference_dataset": "OneScience/QCPINN-Cavity"
+}

config/config.yaml

@@ -0,0 +1,62 @@
+experiment: QCPINN_Table3_Cavity_Angle_Cascade
+paper: https://arxiv.org/pdf/2503.16678
+paths:
+  reference: data/cavity.mat
+  results: results
+  weights: weight
+model:
+  input_dim: 3
+  output_dim: 3
+  hidden_dim: 50
+  n_qubits: 5
+  quantum_layers: 1
+  expected_parameters: 923
+  dtype: float32
+  device: cpu
+  backend: default.qubit
+  shots: null
+physics:
+  rho: 1056.0
+  mu: 0.01
+  t_max: 10.0
+  loss_formula: paper_literal
+  loss_weights: {physics: 0.1, lid: 2.0, walls: 2.0, initial: 4.0}
+training:
+  epochs: 20000
+  batch_size: 64
+  seeds: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
+  optimizer: Adam
+  learning_rate: 0.005
+  betas: [0.9, 0.999]
+  epsilon: 1.0e-8
+  weight_decay: 0.0
+  scheduler: {name: ReduceLROnPlateau, factor: 0.9, patience: 1000, threshold: 0.0001}
+  gradient_clip: 1.0
+  early_stopping: false
+  num_threads: 1
+  checkpoint_every: 1000
+sampling:
+  method: Sobol
+  scramble: true
+  interior: 16384
+  per_wall: 4096
+  initial: 4096
+validation:
+  reference_key: cavity_internal
+  size: 1024
+  seed: 20260325
+  every: 10
+  batch_size: 512
+  selection: mean_relative_l2_u_v_p
+inference:
+  batch_size: 4096
+  reference_keys: [cavity_internal, cavity_gamma0, cavity_gamma1]
+visualization:
+  times: [1.0, 5.0, 10.0]
+  dpi: 200
+assumptions:
+  - Sobol pool sizes, scrambling, seeds 0..9 and validation split are not specified by the paper.
+  - Use literal printed Appendix A.2 losses; component cancellation is a known limitation.
+  - Use independent CRX parameters per Table 2; do not reproduce official parameter reuse bug.
+  - Input coordinates and reference pressure are not normalized, shifted or rescaled.
+  - One epoch is one independently sampled batch of 64 per constraint and one Adam update.

model/QCPINN.py

@@ -0,0 +1,53 @@
+"""Independent implementation of arXiv:2503.16678, Fig.2b and Tables 2/3.
+
+Inputs are physical (t,x,y), outputs (u,v,p). No reference fields enter forward.
+The three groups of n circuit angles are independent, as specified by 3*n*L.
+"""
+import pennylane as qml
+import torch
+from torch import nn
+
+
+class QCPINN(nn.Module):
+    def __init__(self, input_dim=3, output_dim=3, hidden_dim=50,
+                 n_qubits=5, quantum_layers=1, backend="default.qubit", shots=None):
+        super().__init__()
+        if shots is not None:
+            raise ValueError("The paper requires analytic simulation: shots=None")
+        self.n_qubits = n_qubits
+        self.quantum_layers = quantum_layers
+        self.preprocessor = nn.Sequential(nn.Linear(input_dim, hidden_dim), nn.Tanh(),
+                                         nn.Linear(hidden_dim, n_qubits))
+        self.quantum_angles = nn.Parameter(torch.empty(quantum_layers, 3*n_qubits))
+        self.postprocessor = nn.Sequential(nn.Linear(n_qubits, hidden_dim), nn.Tanh(),
+                                          nn.Linear(hidden_dim, output_dim))
+        for block in (self.preprocessor, self.postprocessor):
+            for module in block:
+                if isinstance(module, nn.Linear):
+                    nn.init.xavier_normal_(module.weight)
+                    nn.init.zeros_(module.bias)
+        nn.init.xavier_normal_(self.quantum_angles)
+        device = qml.device(backend, wires=n_qubits, shots=None)
+
+        @qml.qnode(device, interface="torch", diff_method="backprop")
+        def circuit(features, angles):
+            qml.AngleEmbedding(features, wires=range(n_qubits), rotation="X")
+            for layer in range(quantum_layers):
+                for wire in range(n_qubits):
+                    qml.RX(angles[layer, wire], wires=wire)
+                for wire in range(n_qubits):
+                    qml.RZ(angles[layer, n_qubits+wire], wires=wire)
+                # Ring order matches Fig.2b, with distinct angles for every CRX.
+                edges = [(n_qubits-1, 0)] + [(i-1, i) for i in range(n_qubits-1, 0, -1)]
+                for index, edge in enumerate(edges):
+                    qml.CRX(angles[layer, 2*n_qubits+index], wires=edge)
+            return tuple(qml.expval(qml.PauliZ(w)) for w in range(n_qubits))
+
+        self.circuit = circuit
+
+    def forward(self, coordinates):
+        if coordinates.ndim != 2 or coordinates.shape[1] != 3:
+            raise ValueError(f"Expected (batch,3) (t,x,y); got {coordinates.shape}")
+        encoded = self.preprocessor(coordinates)
+        features = torch.stack(self.circuit(encoded, self.quantum_angles), dim=-1)
+        return self.postprocessor(features.to(dtype=coordinates.dtype))

model/__init__.py

@@ -0,0 +1 @@
+from .QCPINN import QCPINN

requirements.txt

@@ -0,0 +1,10 @@
+torch==2.6.0+cpu
+pennylane==0.40.0
+autoray==0.6.12
+numpy==1.26.4
+scipy==1.15.2
+h5py
+matplotlib
+PyYAML
+jax==0.4.35
+jaxlib==0.4.35

results/best_checkpoint_selection.json

@@ -0,0 +1,30 @@
+{
+  "seed": 6,
+  "epoch": 1,
+  "selection_score": 1.3976091685792522,
+  "metrics": {
+    "u": {
+      "relative_l2": 1.2418883137227157,
+      "relative_l2_percent": 124.18883137227157,
+      "error_sum_squares": 67.4534003712654,
+      "reference_sum_squares": 43.73597039713087,
+      "reason": null
+    },
+    "v": {
+      "relative_l2": 1.2429718585473652,
+      "relative_l2_percent": 124.29718585473653,
+      "error_sum_squares": 27.36086485363873,
+      "reference_sum_squares": 17.709537880486465,
+      "reason": null
+    },
+    "p": {
+      "relative_l2": 1.7079673334676762,
+      "relative_l2_percent": 170.7967333467676,
+      "error_sum_squares": 22.096531139899867,
+      "reference_sum_squares": 7.574692034445659,
+      "reason": null
+    }
+  },
+  "checkpoint": "weight/best_model.pt",
+  "phase": "all_ten_formal_runs_completed"
+}

results/best_model/metrics.json

@@ -0,0 +1,168 @@
+{
+  "checkpoint": "weight/best_model.pt",
+  "checkpoint_sha256": "f24e5f3357b045ed2a49fdb8a377cd3668a5cbb8100ff3c13906e3bfce9aff11",
+  "checkpoint_sha256_scope": "original_training_file; published_copy_changes_reference_path_metadata_only",
+  "source": "data/cavity.mat",
+  "source_sha256": "4010d55dc23aa0f8ba79b24baf75e3903044bc4333d772168121cd30bc95c3ba",
+  "seed": 6,
+  "epoch": 1,
+  "selection_score": 1.3976091685792522,
+  "validation_metrics": {
+    "u": {
+      "relative_l2": 1.2418883137227157,
+      "relative_l2_percent": 124.18883137227157,
+      "error_sum_squares": 67.4534003712654,
+      "reference_sum_squares": 43.73597039713087,
+      "reason": null
+    },
+    "v": {
+      "relative_l2": 1.2429718585473652,
+      "relative_l2_percent": 124.29718585473653,
+      "error_sum_squares": 27.36086485363873,
+      "reference_sum_squares": 17.709537880486465,
+      "reason": null
+    },
+    "p": {
+      "relative_l2": 1.7079673334676762,
+      "relative_l2_percent": 170.7967333467676,
+      "error_sum_squares": 22.096531139899867,
+      "reference_sum_squares": 7.574692034445659,
+      "reason": null
+    }
+  },
+  "metric": "raw ||pred-ref||2/||ref||2; no pressure gauge alignment",
+  "groups": {
+    "cavity_internal": {
+      "count": 1010000,
+      "metrics": {
+        "u": {
+          "relative_l2": 1.2369659153068235,
+          "relative_l2_percent": 123.69659153068235,
+          "error_sum_squares": 64793.906059837216,
+          "reference_sum_squares": 42346.61459707969,
+          "reason": null
+        },
+        "v": {
+          "relative_l2": 1.2143086629586421,
+          "relative_l2_percent": 121.43086629586422,
+          "error_sum_squares": 27837.676390170618,
+          "reference_sum_squares": 18878.817807850282,
+          "reason": null
+        },
+        "p": {
+          "relative_l2": 1.6418813661092047,
+          "relative_l2_percent": 164.18813661092045,
+          "error_sum_squares": 20744.156324328287,
+          "reference_sum_squares": 7695.063862736005,
+          "reason": null
+        }
+      }
+    },
+    "cavity_gamma0": {
+      "count": 30300,
+      "metrics": {
+        "u": {
+          "relative_l2": null,
+          "relative_l2_percent": null,
+          "error_sum_squares": 721.6388498068019,
+          "reference_sum_squares": 0.0,
+          "reason": "zero_reference_norm"
+        },
+        "v": {
+          "relative_l2": null,
+          "relative_l2_percent": null,
+          "error_sum_squares": 299.3688306501471,
+          "reference_sum_squares": 0.0,
+          "reason": "zero_reference_norm"
+        },
+        "p": {
+          "relative_l2": 1.0554335687908523,
+          "relative_l2_percent": 105.54335687908522,
+          "error_sum_squares": 2809.0431303861237,
+          "reference_sum_squares": 2521.718480946791,
+          "reason": null
+        }
+      }
+    },
+    "cavity_gamma1": {
+      "count": 10100,
+      "metrics": {
+        "u": {
+          "relative_l2": 0.8486121910397583,
+          "relative_l2_percent": 84.86121910397583,
+          "error_sum_squares": 7273.440772891121,
+          "reference_sum_squares": 10100.0,
+          "reason": null
+        },
+        "v": {
+          "relative_l2": null,
+          "relative_l2_percent": null,
+          "error_sum_squares": 75.20828866056323,
+          "reference_sum_squares": 0.0,
+          "reason": "zero_reference_norm"
+        },
+        "p": {
+          "relative_l2": 1.003925357845928,
+          "relative_l2_percent": 100.3925357845928,
+          "error_sum_squares": 2125.860076265305,
+          "reference_sum_squares": 2109.2683099242454,
+          "reason": null
+        }
+      }
+    }
+  },
+  "reference_regenerated": false,
+  "checkpoint_reloaded": true,
+  "internal_holdout": {
+    "count": 1008976,
+    "metrics": {
+      "u": {
+        "relative_l2": 1.2369608160205556,
+        "relative_l2_percent": 123.69608160205556,
+        "error_sum_squares": 64726.45265992043,
+        "reference_sum_squares": 42302.87862668262,
+        "reason": null
+      },
+      "v": {
+        "relative_l2": 1.2142814318609751,
+        "relative_l2_percent": 121.42814318609751,
+        "error_sum_squares": 27810.31552531313,
+        "reference_sum_squares": 18861.10826996979,
+        "reason": null
+      },
+      "p": {
+        "relative_l2": 1.6418149379825397,
+        "relative_l2_percent": 164.18149379825397,
+        "error_sum_squares": 20722.05979287031,
+        "reference_sum_squares": 7687.489170701568,
+        "reason": null
+      }
+    }
+  },
+  "all_reference": {
+    "metrics": {
+      "u": {
+        "relative_l2": 1.1780781571826924,
+        "relative_l2_percent": 117.80781571826924,
+        "error_sum_squares": 72788.98568253515,
+        "reference_sum_squares": 52446.61459707969,
+        "reason": null
+      },
+      "v": {
+        "relative_l2": 1.222451087407963,
+        "relative_l2_percent": 122.2451087407963,
+        "error_sum_squares": 28212.25350948133,
+        "reference_sum_squares": 18878.817807850282,
+        "reason": null
+      },
+      "p": {
+        "relative_l2": 1.44336969627177,
+        "relative_l2_percent": 144.33696962717698,
+        "error_sum_squares": 25679.059530979714,
+        "reference_sum_squares": 12326.050653607042,
+        "reason": null
+      }
+    },
+    "count": 1050400
+  }
+}

results/best_model/report.md

@@ -0,0 +1,43 @@
+执行结论：完整运行已完成，但未复现Table 3数值精度。全部训练仅使用本案例新提交的计算作业975143；没有借用、取消或修改账号中的既有作业。10个种子各完成20000epoch，正式权重重载与可视化均已完成。
+
+按论文印刷式逐字实现的loss允许u=.5,v=-.5,p=0的零损失非物理解。本实验保留这一公式，未用官方额外监督项替换它，也未另造参考真值；来源冲突与实现假设详见[审计规格](../../reproduction_spec.md)。
+
+# QCPINN Angle-Cascade cavity reproduction
+
+Selected seed: 6; epoch: 1; 923 parameters.
+
+Reference values were read directly from the supplied cavity.mat; no generated reference truth.
+
+| Field | Internal relative L2 (%) | Paper Table 3 mean ± std (%) |
+|---|---:|---:|
+| u | 123.696592 | 25.73 ± 2.49 |
+| v | 121.430866 | 38.00 ± 3.76 |
+| p | 164.188137 | 42.65 ± 8.69 |
+
+## Ten independent complete runs (20000 epochs each)
+
+| Field | Final checkpoints mean ± std (%) | Validation-best checkpoints mean ± std (%) | Paper Table 3 (%) |
+|---|---:|---:|---:|
+| u | 263.827208 ± 0.072903 | 200.415311 ± 61.786860 | 25.73 ± 2.49 |
+| v | 379.130934 ± 0.111758 | 274.729384 ± 87.551613 | 38.00 ± 3.76 |
+| p | 100.122990 ± 0.609755 | 148.862860 ± 73.002451 | 42.65 ± 8.69 |
+
+Final loss: 4.71713587e-07 ± 4.71234223e-07; paper: 0.10 ± 0.02.
+Time per update: 0.216483 ± 0.003915 s; wall time for PDE/BC/IC forward, backward, gradient clip and Adam step; validation/checkpoint I/O excluded.
+Peak process RSS: 712.232 ± 28.112 MB. Actual hardware is recorded in results/compute_hardware.txt; timing/memory are not hardware-identical to the paper.
+Population standard deviation (ddof=0) is used; the paper does not specify this convention.
+
+The selected best checkpoint is a single model; the paper entries are ten-run statistics.
+See table3_comparison.json for separately aggregated final and best checkpoints.
+
+## Explicit limitations
+- Sobol pool sizes, scrambling, seeds 0..9 and validation split are not specified by the paper.
+- Use literal printed Appendix A.2 losses; component cancellation is a known limitation.
+- Use independent CRX parameters per Table 2; do not reproduce official parameter reuse bug.
+- Input coordinates and reference pressure are not normalized, shifted or rescaled.
+- One epoch is one independently sampled batch of 64 per constraint and one Adam update.
+- Appendix A.2 printed sums can cancel between components. The literal printed losses were retained.
+- The reference contains saved time intervals of 0.1 s, while the paper describes a solver step of 0.01 s.
+- Validation uses a fixed reference subset for checkpoint selection. Holdout metrics exclude those records.
+
+Source: [paper](https://arxiv.org/pdf/2503.16678), sections 4/5, Figure 2, Tables 2/3 and Appendix A.2.

results/compute_hardware.txt

@@ -0,0 +1,25 @@
+a01r3n10
+Architecture:        x86_64
+CPU op-mode(s):      32-bit, 64-bit
+Byte Order:          Little Endian
+CPU(s):              128
+On-line CPU(s) list: 0-127
+Thread(s) per core:  1
+Core(s) per socket:  64
+Socket(s):           2
+NUMA node(s):        8
+Vendor ID:           HygonGenuine
+CPU family:          24
+Model:               7
+Model name:          Hygon C86 Processor
+Stepping:            0
+CPU MHz:             2400.017
+BogoMIPS:            4800.03
+Virtualization:      AMD-V
+L1d cache:           32K
+L1i cache:           32K
+L2 cache:            512K
+L3 cache:            16384K
+NUMA node0 CPU(s):   0-15
+NUMA node1 CPU(s):   16-31
+NUMA node2 CPU(s):   32-47

results/table3_comparison.json

@@ -0,0 +1,179 @@
+{
+  "status": "completed",
+  "runs": 10,
+  "epochs_per_run": 20000,
+  "statistics": {
+    "final": {
+      "u": {
+        "mean_percent": 263.82720846357546,
+        "std_percent": 0.07290313116389122,
+        "values_percent": [
+          263.8550091118908,
+          263.9313754674022,
+          263.8202844356478,
+          263.7872862846581,
+          263.86905668107204,
+          263.7875723279372,
+          263.8015990771984,
+          263.9572048027661,
+          263.7136317704975,
+          263.74906467668444
+        ]
+      },
+      "v": {
+        "mean_percent": 379.1309335658924,
+        "std_percent": 0.11175777875192125,
+        "values_percent": [
+          379.0429768655644,
+          379.2406101730609,
+          379.0715007982627,
+          379.04294338051415,
+          379.1214781463276,
+          379.04151966547096,
+          379.3142254787664,
+          379.2905145628615,
+          379.17106663969815,
+          378.9724999483974
+        ]
+      },
+      "p": {
+        "mean_percent": 100.12298999068985,
+        "std_percent": 0.6097553354196278,
+        "values_percent": [
+          99.99971594555477,
+          99.96530609848863,
+          101.92920783109811,
+          99.79079939688157,
+          99.76182788947703,
+          99.81738607251161,
+          100.08049441555151,
+          99.93770620014934,
+          99.95562455821036,
+          99.9918314989756
+        ]
+      }
+    },
+    "best": {
+      "u": {
+        "mean_percent": 200.41531063983444,
+        "std_percent": 61.786859912537416,
+        "values_percent": [
+          216.73178735741482,
+          257.8566730563977,
+          244.05346210075695,
+          223.0109183184541,
+          209.36992710804304,
+          98.65206632622584,
+          123.69659153068235,
+          261.36667138461894,
+          108.33561949271215,
+          261.0793897230385
+        ]
+      },
+      "v": {
+        "mean_percent": 274.72938385999566,
+        "std_percent": 87.55161329864875,
+        "values_percent": [
+          306.2224081242297,
+          348.16036507422353,
+          351.4166323499106,
+          308.4572411231171,
+          299.9983810702734,
+          166.83810214589104,
+          121.43086629586422,
+          339.46700708938135,
+          147.04576328216456,
+          358.25707204490146
+        ]
+      },
+      "p": {
+        "mean_percent": 148.86286014952978,
+        "std_percent": 73.00245075399947,
+        "values_percent": [
+          103.59704594765071,
+          100.38273398618911,
+          99.33596465585998,
+          113.92679231218537,
+          128.94455232416792,
+          227.51965246296444,
+          164.18813661092045,
+          111.02781481911512,
+          335.8090833036171,
+          103.89682507262745
+        ]
+      }
+    }
+  },
+  "std_ddof": 0,
+  "paper_table3_percent": {
+    "u": {
+      "mean": 25.73,
+      "std": 2.49
+    },
+    "v": {
+      "mean": 38.0,
+      "std": 3.76
+    },
+    "p": {
+      "mean": 42.65,
+      "std": 8.69
+    }
+  },
+  "final_loss": {
+    "mean": 4.7171358659170436e-07,
+    "std": 4.712342234329683e-07
+  },
+  "parameters": 923,
+  "time_per_update_seconds": {
+    "mean": 0.21648326783038208,
+    "std": 0.003914749140125143,
+    "run_means": [
+      0.21697817574576475,
+      0.21848320898595266,
+      0.2186092604946112,
+      0.22009400647231378,
+      0.22067649887888693,
+      0.2170402910890989,
+      0.21741632961235008,
+      0.2171733768468257,
+      0.21084768179748206,
+      0.2075138483805349
+    ],
+    "definition": "wall time for PDE/BC/IC forward, backward, gradient clip and Adam step; validation/checkpoint I/O excluded"
+  },
+  "peak_process_rss_mb": {
+    "mean": 712.231640625,
+    "std": 28.112347417123058,
+    "values": [
+      676.7890625,
+      680.02734375,
+      677.796875,
+      676.68359375,
+      735.2109375,
+      734.66796875,
+      735.0390625,
+      736.55859375,
+      735.4140625,
+      734.12890625
+    ]
+  },
+  "paper_table3_other": {
+    "parameters": 923,
+    "loss_mean": 0.1,
+    "loss_std": 0.02,
+    "time_seconds_mean": 0.36,
+    "time_seconds_std": 0.01,
+    "memory_mb_mean": 540.0,
+    "memory_mb_std": 0.0
+  },
+  "selected_seed": 6,
+  "selection_metric": "mean validation relative L2",
+  "best_model": "weight/best_model.pt",
+  "limitations": [
+    "Sobol pool sizes, scrambling, seeds 0..9 and validation split are not specified by the paper.",
+    "Use literal printed Appendix A.2 losses; component cancellation is a known limitation.",
+    "Use independent CRX parameters per Table 2; do not reproduce official parameter reuse bug.",
+    "Input coordinates and reference pressure are not normalized, shifted or rescaled.",
+    "One epoch is one independently sampled batch of 64 per constraint and one Adam update."
+  ]
+}

scripts/common.py

@@ -0,0 +1,198 @@
+"""Data, Sobol, PDE and metric contracts shared by the three entry points."""
+import hashlib
+import json
+import os
+import random
+import socket
+from pathlib import Path
+
+import h5py
+import numpy as np
+import torch
+import yaml
+
+ROOT = Path(__file__).resolve().parents[1]
+
+
+def load_config(path):
+    with open(path, encoding="utf-8") as stream:
+        config = yaml.safe_load(stream)
+    if config["training"]["batch_size"] != 64 or config["training"]["epochs"] != 20000:
+        raise ValueError("Formal protocol must retain batch_size=64 and epochs=20000")
+    if config["training"]["early_stopping"]:
+        raise ValueError("User requires complete formal training; early stopping disabled")
+    return config
+
+
+def resolved(config, name):
+    path = Path(config["paths"][name])
+    return path if path.is_absolute() else ROOT / path
+
+
+def atomic_json(path, payload):
+    path = Path(path)
+    path.parent.mkdir(parents=True, exist_ok=True)
+    tmp = path.with_suffix(path.suffix + ".tmp")
+    tmp.write_text(json.dumps(payload, ensure_ascii=False, indent=2, allow_nan=False), encoding="utf-8")
+    tmp.replace(path)
+
+
+def atomic_save(path, payload):
+    path = Path(path)
+    path.parent.mkdir(parents=True, exist_ok=True)
+    tmp = path.with_suffix(path.suffix + ".tmp")
+    torch.save(payload, tmp)
+    tmp.replace(path)
+
+
+def require_compute_node():
+    # A self-submitted allocation also records a job name unique to this case.
+    jid, name = os.environ.get("SLURM_JOB_ID"), os.environ.get("SLURM_JOB_NAME", "")
+    if not jid or not name.startswith("qcpinn_cavity_") or "login" in socket.gethostname().lower():
+        raise RuntimeError("Training/smoke requires a new case-owned SLURM compute allocation")
+    return {"job_id": jid, "job_name": name, "hostname": socket.gethostname(),
+            "array_job_id": os.environ.get("SLURM_ARRAY_JOB_ID"),
+            "array_task_id": os.environ.get("SLURM_ARRAY_TASK_ID")}
+
+
+def seed_everything(seed, threads=1):
+    random.seed(seed)
+    np.random.seed(seed)
+    torch.manual_seed(seed)
+    torch.set_num_threads(threads)
+
+
+def build_model(config):
+    from model.QCPINN import QCPINN
+    m = config["model"]
+    model = QCPINN(**{key: m[key] for key in ["input_dim", "output_dim", "hidden_dim", "n_qubits", "quantum_layers", "backend", "shots"]})
+    model = model.to(dtype=getattr(torch, m["dtype"]), device=m["device"])
+    n = sum(p.numel() for p in model.parameters())
+    if n != m["expected_parameters"]:
+        raise ValueError(f"Parameter count {n} != paper count {m['expected_parameters']}")
+    return model
+
+
+def make_sobol_pools(config, seed):
+    s = config["sampling"]
+    pools = {}
+    for index, kind in enumerate(["interior", "lid", "left", "right", "bottom", "initial"]):
+        count = s["interior"] if kind == "interior" else s["initial"] if kind == "initial" else s["per_wall"]
+        engine = torch.quasirandom.SobolEngine(3, scramble=s["scramble"], seed=seed*17+index)
+        xyz = engine.draw(count).to(dtype=getattr(torch, config["model"]["dtype"]))
+        xyz[:, 0] *= config["physics"]["t_max"]
+        if kind == "lid": xyz[:, 2] = 1.0
+        if kind == "left": xyz[:, 1] = 0.0
+        if kind == "right": xyz[:, 1] = 1.0
+        if kind == "bottom": xyz[:, 2] = 0.0
+        if kind == "initial": xyz[:, 0] = 0.0
+        pools[kind] = xyz.to(config["model"]["device"])
+    return pools
+
+
+def sample_batches(pools, size, rng):
+    return {key: points[torch.as_tensor(rng.choice(len(points), size=size, replace=False), device=points.device)]
+            for key, points in pools.items()}
+
+
+def derivative(field, coords):
+    return torch.autograd.grad(field, coords, torch.ones_like(field), create_graph=True)[0]
+
+
+def pde_residuals(model, coordinates, physics):
+    coords = coordinates.detach().clone().requires_grad_(True)
+    uvp = model(coords)
+    u, v, p = uvp.unbind(-1)
+    du, dv, dp = [derivative(field, coords) for field in (u, v, p)]
+    uxx = derivative(du[:, 1], coords)[:, 1]
+    uyy = derivative(du[:, 2], coords)[:, 2]
+    vxx = derivative(dv[:, 1], coords)[:, 1]
+    vyy = derivative(dv[:, 2], coords)[:, 2]
+    ru = du[:, 0] + u*du[:, 1] + v*du[:, 2] + dp[:, 1]/physics["rho"] - physics["mu"]*(uxx+uyy)
+    rv = dv[:, 0] + u*dv[:, 1] + v*dv[:, 2] + dp[:, 2]/physics["rho"] - physics["mu"]*(vxx+vyy)
+    rc = du[:, 1] + dv[:, 2]
+    return torch.stack((ru, rv, rc), dim=-1)
+
+
+def loss_terms(model, batches, config):
+    physics = config["physics"]
+    if physics["loss_formula"] != "paper_literal":
+        raise ValueError("This reproduction implements the literal Appendix A.2 losses")
+    residuals = pde_residuals(model, batches["interior"], physics)
+    terms = {"physics": residuals.square().mean(0).sum()}
+    names = ["lid", "left", "right", "bottom", "initial"]
+    sizes = [len(batches[k]) for k in names]
+    fields = torch.split(model(torch.cat([batches[k] for k in names])), sizes)
+    lid, *rest = fields
+    terms["lid"] = (1-lid[:, 0]+lid[:, 1]).square().mean()
+    terms["walls"] = sum(out[:, :2].sum(-1).square().mean() for out in rest[:3])
+    terms["initial"] = rest[3].sum(-1).square().mean()
+    terms["total"] = sum(physics["loss_weights"][key]*value for key, value in terms.items())
+    return terms
+
+
+def validation_indices(config, total):
+    size = config["validation"]["size"]
+    if size >= total:
+        raise ValueError("Validation must leave nonempty holdout")
+    return np.sort(np.random.default_rng(config["validation"]["seed"]).choice(total, size, replace=False))
+
+
+def load_validation(config):
+    key = config["validation"]["reference_key"]
+    with h5py.File(resolved(config, "reference"), "r") as source:
+        data = source[key]
+        if data.shape[0] != 6:
+            raise ValueError("Expected MATLAB v7.3 dataset of shape (6,N)")
+        ids = validation_indices(config, data.shape[1])
+        rows = data[:, ids].T
+    if not np.isfinite(rows).all():
+        raise ValueError("Reference validation contains NaN/Inf")
+    return rows[:, :3], rows[:, 3:], ids
+
+
+def relative_metrics(error_squares, reference_squares):
+    metrics = {}
+    for i, key in enumerate(("u", "v", "p")):
+        den, num = float(reference_squares[i]), float(error_squares[i])
+        value = float(np.sqrt(num/den)) if den > 0 else None
+        metrics[key] = {"relative_l2": value, "relative_l2_percent": None if value is None else 100*value,
+                        "error_sum_squares": num, "reference_sum_squares": den,
+                        "reason": None if value is not None else "zero_reference_norm"}
+    return metrics
+
+
+@torch.no_grad()
+def predict(model, coords, size):
+    param = next(model.parameters())
+    return np.concatenate([model(torch.as_tensor(coords[i:i+size], dtype=param.dtype, device=param.device)).cpu().numpy()
+                           for i in range(0, len(coords), size)], axis=0)
+
+
+def evaluate(model, coordinates, reference, size):
+    was_training = model.training
+    model.eval()
+    predicted = predict(model, coordinates, size).astype(np.float64)
+    model.train(was_training)
+    errors = predicted-reference
+    return relative_metrics(np.sum(errors**2, 0), np.sum(reference**2, 0))
+
+
+def selection_score(metrics):
+    values = [metrics[k]["relative_l2"] for k in ("u", "v", "p")]
+    if any(v is None or not np.isfinite(v) for v in values):
+        raise ValueError("Validation L2 is undefined/nonfinite")
+    return float(np.mean(values))
+
+
+def sha256(path):
+    digest = hashlib.sha256()
+    with open(path, "rb") as stream:
+        for block in iter(lambda: stream.read(1024*1024), b""): digest.update(block)
+    return digest.hexdigest()
+
+
+def load_checkpoint(path, model):
+    checkpoint = torch.load(path, map_location=next(model.parameters()).device, weights_only=False)
+    model.load_state_dict(checkpoint["model_state_dict"], strict=True)
+    return checkpoint

scripts/inference.py

@@ -0,0 +1,98 @@
+"""Reload a trained checkpoint; stream evaluation from the supplied MATLAB HDF5."""
+import argparse
+import json
+import sys
+from pathlib import Path
+sys.path.insert(0, str(Path(__file__).resolve().parents[1]))
+import h5py
+import numpy as np
+import torch
+from common import (ROOT, atomic_json, build_model, load_checkpoint, load_config,
+                    predict, relative_metrics, resolved, seed_everything, sha256,
+                    validation_indices)
+
+
+def infer(config, checkpoint_path, output_dir, metrics_only=False):
+    seed_everything(0, config["training"]["num_threads"])
+    model = build_model(config)
+    checkpoint = load_checkpoint(checkpoint_path, model)
+    if checkpoint["config"]["model"] != config["model"]:
+        raise ValueError("Checkpoint and requested model configuration differ")
+    model.eval()
+    out = Path(output_dir)
+    out.mkdir(parents=True, exist_ok=True)
+    source_path = resolved(config, "reference")
+    report = {"checkpoint": str(Path(checkpoint_path).resolve()), "checkpoint_sha256": sha256(checkpoint_path),
+              "source": str(source_path), "source_sha256": sha256(source_path),
+              "seed": checkpoint["seed"], "epoch": checkpoint["epoch"], "selection_score": checkpoint["selection_score"],
+              "validation_metrics": checkpoint["validation_metrics"],
+              "metric": "raw ||pred-ref||2/||ref||2; no pressure gauge alignment",
+              "groups": {}, "reference_regenerated": False, "checkpoint_reloaded": True}
+    size = config["inference"]["batch_size"]
+    global_num, global_den = np.zeros(3), np.zeros(3)
+    output = None if metrics_only else h5py.File(out/"predictions.h5.tmp", "w")
+    try:
+        with h5py.File(source_path, "r") as source:
+            for key in config["inference"]["reference_keys"]:
+                data = source[key]
+                if data.ndim != 2 or data.shape[0] != 6: raise ValueError(f"Invalid {key} shape: {data.shape}")
+                count = data.shape[1]
+                numer, denom, hold_num, hold_den = [np.zeros(3) for _ in range(4)]
+                if key == config["validation"]["reference_key"]:
+                    valids = validation_indices(config, count)
+                group = None
+                if output is not None:
+                    group = output.create_group(key)
+                    group.create_dataset("coordinates", (count, 3), dtype="f8", chunks=True)
+                    group.create_dataset("prediction", (count, 3), dtype="f4", chunks=True)
+                for start in range(0, count, size):
+                    end = min(start+size, count)
+                    rows = data[:, start:end].T
+                    if not np.isfinite(rows).all(): raise ValueError(f"Nonfinite input in {key} at {start}")
+                    preds = predict(model, rows[:, :3], size)
+                    if not np.isfinite(preds).all(): raise FloatingPointError("Nonfinite predictions")
+                    error2 = (preds.astype(np.float64)-rows[:, 3:])**2
+                    ref2 = rows[:, 3:]**2
+                    numer += error2.sum(0)
+                    denom += ref2.sum(0)
+                    if key == config["validation"]["reference_key"]:
+                        keep = ~np.isin(np.arange(start, end), valids, assume_unique=True)
+                        hold_num += error2[keep].sum(0)
+                        hold_den += ref2[keep].sum(0)
+                    if group is not None:
+                        group["coordinates"][start:end] = rows[:, :3]
+                        group["prediction"][start:end] = preds
+                metrics = relative_metrics(numer, denom)
+                report["groups"][key] = {"count": count, "metrics": metrics}
+                global_num += numer
+                global_den += denom
+                if key == config["validation"]["reference_key"]:
+                    report["internal_holdout"] = {"count": count-len(valids), "metrics": relative_metrics(hold_num, hold_den)}
+                print(f"INFERENCE key={key} count={count} metrics={json.dumps(metrics)}", flush=True)
+        report["all_reference"] = {"metrics": relative_metrics(global_num, global_den),
+                                    "count": sum(v["count"] for v in report["groups"].values())}
+        if output is not None:
+            output.attrs["provenance"] = json.dumps({k: v for k, v in report.items() if k not in ["groups", "internal_holdout", "all_reference"]})
+            output.close()
+            output = None
+            (out/"predictions.h5.tmp").replace(out/"predictions.h5")
+        atomic_json(out/"metrics.json", report)
+    finally:
+        if output is not None: output.close()
+    print("INFERENCE_COMPLETED", flush=True)
+    return report
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("--config", default=str(ROOT/"config/config.yaml"))
+    parser.add_argument("--checkpoint", type=Path)
+    parser.add_argument("--output-dir", type=Path)
+    parser.add_argument("--seed", type=int)
+    parser.add_argument("--kind", choices=["best", "last"], default="best")
+    parser.add_argument("--metrics-only", action="store_true")
+    args = parser.parse_args()
+    cfg = load_config(args.config)
+    checkpoint = args.checkpoint or (resolved(cfg, "weights")/f"seed_{args.seed}"/("best_model.pt" if args.kind == "best" else "last_model.pt") if args.seed is not None else resolved(cfg, "weights")/"best_model.pt")
+    output = args.output_dir or (resolved(cfg, "results")/f"seed_{args.seed}"/("best_eval" if args.kind == "best" else "final_eval") if args.seed is not None else resolved(cfg, "results")/"best_model")
+    infer(cfg, checkpoint, output, args.metrics_only)

scripts/result.py

@@ -0,0 +1,188 @@
+"""Create scientific figures and aggregate ten independent Table 3 runs."""
+import argparse
+import csv
+import json
+import shutil
+import sys
+from pathlib import Path
+sys.path.insert(0, str(Path(__file__).resolve().parents[1]))
+import h5py
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import numpy as np
+from common import ROOT, atomic_json, load_config, resolved
+
+
+def aggregate(config):
+    results, weights = resolved(config, "results"), resolved(config, "weights")
+    records = []
+    for seed in config["training"]["seeds"]:
+        directory = results/f"seed_{seed}"
+        summary = json.loads((directory/"training_summary.json").read_text())
+        if summary["completed_epochs"] != 20000: raise ValueError(f"Incomplete seed {seed}")
+        best = json.loads((directory/"best_eval/metrics.json").read_text())
+        final = json.loads((directory/"final_eval/metrics.json").read_text())
+        records.append({"seed": seed, "best": best, "final": final, "training": summary})
+    winner = min(records, key=lambda row: row["best"]["selection_score"])
+    shutil.copy2(weights/f"seed_{winner['seed']}"/"best_model.pt", weights/"best_model.pt")
+    atomic_json(results/"best_checkpoint_selection.json", {"seed": winner["seed"],
+                "epoch": winner["best"]["epoch"], "selection_score": winner["best"]["selection_score"],
+                "metrics": winner["best"]["validation_metrics"],
+                "checkpoint": str(weights/"best_model.pt"), "phase": "all_ten_formal_runs_completed"})
+    reference = {"u": {"mean": 25.73, "std": 2.49}, "v": {"mean": 38.00, "std": 3.76}, "p": {"mean": 42.65, "std": 8.69}}
+    statistics = {}
+    for kind in ["final", "best"]:
+        statistics[kind] = {}
+        for field in ("u", "v", "p"):
+            values = np.array([r[kind]["groups"]["cavity_internal"]["metrics"][field]["relative_l2_percent"] for r in records])
+            statistics[kind][field] = {"mean_percent": float(values.mean()), "std_percent": float(values.std(ddof=0)), "values_percent": values.tolist()}
+    losses = np.array([r["training"]["last_loss"]["total"] for r in records])
+    run_step_times = []
+    for row in records:
+        with open(results/f"seed_{row['seed']}"/"history.csv") as stream:
+            times = [float(item["step_seconds"]) for item in csv.DictReader(stream)]
+        if len(times) != 20000:
+            raise ValueError(f"Expected exactly 20000 recorded updates for seed {row['seed']}")
+        run_step_times.append(float(np.mean(times)))
+    peak_rss = np.array([r["training"]["peak_rss_mb"] for r in records])
+    report = {"status": "completed", "runs": len(records), "epochs_per_run": 20000,
+              "statistics": statistics, "std_ddof": 0, "paper_table3_percent": reference,
+              "final_loss": {"mean": float(losses.mean()), "std": float(losses.std())},
+              "parameters": 923,
+              "time_per_update_seconds": {"mean": float(np.mean(run_step_times)), "std": float(np.std(run_step_times)), "run_means": run_step_times, "definition": "wall time for PDE/BC/IC forward, backward, gradient clip and Adam step; validation/checkpoint I/O excluded"},
+              "peak_process_rss_mb": {"mean": float(peak_rss.mean()), "std": float(peak_rss.std()), "values": peak_rss.tolist()},
+              "paper_table3_other": {"parameters": 923, "loss_mean": 0.10, "loss_std": 0.02, "time_seconds_mean": 0.36, "time_seconds_std": 0.01, "memory_mb_mean": 540.0, "memory_mb_std": 0.0},
+              "selected_seed": winner["seed"], "selection_metric": "mean validation relative L2",
+              "best_model": str(weights/"best_model.pt"),
+              "limitations": config["assumptions"]}
+    atomic_json(results/"table3_comparison.json", report)
+    print(json.dumps(report, indent=2), flush=True)
+    return report
+
+
+def plot_history(config, seed, destination):
+    path = resolved(config, "results")/f"seed_{seed}"/"history.csv"
+    with open(path) as stream: rows = list(csv.DictReader(stream))
+    epochs = np.array([int(r["epoch"]) for r in rows])
+    fig, axes = plt.subplots(1, 2, figsize=(12, 4.2), constrained_layout=True)
+    for key in ("total", "physics", "lid", "walls", "initial"):
+        axes[0].semilogy(epochs, [max(float(r[key]), 1e-30) for r in rows], label=key, linewidth=.8)
+    axes[0].set(xlabel="Epoch (one Adam update)", ylabel="Loss", title=f"Angle-Cascade, seed={seed}")
+    valid = [r for r in rows if r["val_u_l2"]]
+    for key in ("u", "v", "p"):
+        axes[1].plot([int(r["epoch"]) for r in valid], [100*float(r[f"val_{key}_l2"]) for r in valid], label=key)
+    axes[1].set(xlabel="Epoch", ylabel="Validation relative L2 (%)", title="Raw u/v/p, no pressure shift")
+    for ax in axes: ax.grid(alpha=.2); ax.legend()
+    target = destination/"training_curves.png"
+    fig.savefig(target, dpi=config["visualization"]["dpi"])
+    plt.close(fig)
+    return str(target)
+
+
+def plot_fields(config, predictions, destination):
+    images = []
+    with h5py.File(resolved(config, "reference"), "r") as source, h5py.File(predictions, "r") as inference:
+        key = "cavity_internal"
+        raw = source[key]
+        checkpoint_meta = json.loads(inference.attrs["provenance"])
+        actual_times = np.unique(raw[0, :])
+        for desired in config["visualization"]["times"]:
+            t = float(actual_times[np.argmin(abs(actual_times-desired))])
+            ids = np.flatnonzero(np.isclose(raw[0, :], t, atol=1e-8, rtol=0))
+            rows = raw[:, ids].T
+            pred = inference[key]["prediction"][ids]
+            xs, ys = np.unique(rows[:, 1]), np.unique(rows[:, 2])
+            if len(xs)*len(ys) != len(rows): raise ValueError("Reference snapshot is not a complete rectangular grid")
+            order = np.lexsort((rows[:, 1], rows[:, 2]))
+            truth = rows[order, 3:].reshape(len(ys), len(xs), 3)
+            predicted = pred[order].reshape(len(ys), len(xs), 3)
+            fig, axes = plt.subplots(3, 3, figsize=(11, 10), constrained_layout=True)
+            for index, field in enumerate(("u", "v", "p")):
+                low, high = min(truth[..., index].min(), predicted[..., index].min()), max(truth[..., index].max(), predicted[..., index].max())
+                if low == high: high = low+1e-12
+                for col, values in enumerate((truth[..., index], predicted[..., index], abs(predicted[..., index]-truth[..., index]))):
+                    cmap = "coolwarm" if index == 2 and col < 2 else "viridis" if col < 2 else "magma"
+                    kwargs = {"vmin": low, "vmax": high} if col < 2 else {"vmin": 0}
+                    artist = axes[index, col].pcolormesh(xs, ys, values, shading="auto", cmap=cmap, **kwargs)
+                    axes[index, col].set(xlabel="x (m)", ylabel="y (m)", title=f"{field}: {['Reference', 'QCPINN', 'Absolute error'][col]}")
+                    axes[index, col].set_aspect("equal")
+                    fig.colorbar(artist, ax=axes[index, col], shrink=.8)
+            fig.suptitle(f"Lid-driven cavity, t={t:g} s | seed={checkpoint_meta['seed']}, epoch={checkpoint_meta['epoch']} | raw fields")
+            target = destination/f"fields_t{t:g}.png"
+            fig.savefig(target, dpi=config["visualization"]["dpi"])
+            plt.close(fig)
+            images.append(str(target))
+            fig, axes = plt.subplots(1, 2, figsize=(10, 4.5), constrained_layout=True)
+            for ax, data, title in zip(axes, (truth, predicted), ("Reference", "QCPINN")):
+                speed = np.linalg.norm(data[..., :2], axis=-1)
+                artist = ax.streamplot(xs, ys, data[..., 0], data[..., 1], color=speed, cmap="viridis", density=1.1)
+                fig.colorbar(artist.lines, ax=ax, label="Speed (m/s)")
+                ax.set(xlabel="x (m)", ylabel="y (m)", title=f"{title}, t={t:g} s", xlim=(0, 1), ylim=(0, 1))
+                ax.set_aspect("equal")
+            target = destination/f"streamlines_t{t:g}.png"
+            fig.savefig(target, dpi=config["visualization"]["dpi"])
+            plt.close(fig)
+            images.append(str(target))
+    return images
+
+
+def make_results(config, seed=None, output_dir=None, predictions=None):
+    results = resolved(config, "results")
+    comparison = None
+    if seed is None:
+        comparison = json.loads((results/"table3_comparison.json").read_text())
+        seed = comparison["selected_seed"]
+        evaluation = results/"best_model"
+    else: evaluation = results/f"seed_{seed}"/"best_eval"
+    output = Path(output_dir) if output_dir else evaluation/"figures"
+    output.mkdir(parents=True, exist_ok=True)
+    paths = [plot_history(config, seed, output)]
+    paths += plot_fields(config, Path(predictions) if predictions else evaluation/"predictions.h5", output)
+    metrics = json.loads((evaluation/"metrics.json").read_text())
+    lines = ["# QCPINN Angle-Cascade cavity reproduction", "", f"Selected seed: {seed}; epoch: {metrics['epoch']}; 923 parameters.",
+             "", "Reference values were read directly from the supplied cavity.mat; no generated reference truth.",
+             "", "| Field | Internal relative L2 (%) | Paper Table 3 mean ± std (%) |", "|---|---:|---:|"]
+    for k, paper in zip(("u", "v", "p"), ("25.73 ± 2.49", "38.00 ± 3.76", "42.65 ± 8.69")):
+        value = metrics["groups"]["cavity_internal"]["metrics"][k]["relative_l2_percent"]
+        lines.append(f"| {k} | {value:.6f} | {paper} |")
+    if comparison is not None:
+        lines += ["", "## Ten independent complete runs (20000 epochs each)", "",
+                  "| Field | Final checkpoints mean ± std (%) | Validation-best checkpoints mean ± std (%) | Paper Table 3 (%) |",
+                  "|---|---:|---:|---:|"]
+        for field in ("u", "v", "p"):
+            final = comparison["statistics"]["final"][field]
+            best = comparison["statistics"]["best"][field]
+            paper = comparison["paper_table3_percent"][field]
+            lines.append(f"| {field} | {final['mean_percent']:.6f} ± {final['std_percent']:.6f} | {best['mean_percent']:.6f} ± {best['std_percent']:.6f} | {paper['mean']:.2f} ± {paper['std']:.2f} |")
+        loss = comparison["final_loss"]
+        timing = comparison["time_per_update_seconds"]
+        memory = comparison["peak_process_rss_mb"]
+        lines += ["", f"Final loss: {loss['mean']:.8e} ± {loss['std']:.8e}; paper: 0.10 ± 0.02.",
+                  f"Time per update: {timing['mean']:.6f} ± {timing['std']:.6f} s; {timing['definition']}.",
+                  f"Peak process RSS: {memory['mean']:.3f} ± {memory['std']:.3f} MB. Actual hardware is recorded in results/compute_hardware.txt; timing/memory are not hardware-identical to the paper.",
+                  "Population standard deviation (ddof=0) is used; the paper does not specify this convention."]
+    lines += ["", "The selected best checkpoint is a single model; the paper entries are ten-run statistics.",
+              "See table3_comparison.json for separately aggregated final and best checkpoints.", "", "## Explicit limitations"]
+    lines += [f"- {item}" for item in config["assumptions"]]
+    lines += ["- Appendix A.2 printed sums can cancel between components. The literal printed losses were retained.",
+              "- The reference contains saved time intervals of 0.1 s, while the paper describes a solver step of 0.01 s.",
+              "- Validation uses a fixed reference subset for checkpoint selection. Holdout metrics exclude those records.",
+              "", "Source: [paper](https://arxiv.org/pdf/2503.16678), sections 4/5, Figure 2, Tables 2/3 and Appendix A.2."]
+    report = evaluation/"report.md"
+    report.write_text("\n".join(lines)+"\n")
+    atomic_json(evaluation/"visualization_manifest.json", {"figures": paths, "report": str(report), "dpi": config["visualization"]["dpi"], "reference": str(resolved(config, "reference"))})
+    print(f"VISUALIZATION_COMPLETED figures={len(paths)} report={report}", flush=True)
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("--config", default=str(ROOT/"config/config.yaml"))
+    parser.add_argument("--aggregate", action="store_true")
+    parser.add_argument("--seed", type=int)
+    parser.add_argument("--output-dir", type=Path)
+    parser.add_argument("--predictions", type=Path)
+    args = parser.parse_args()
+    cfg = load_config(args.config)
+    if args.aggregate: aggregate(cfg)
+    else: make_results(cfg, args.seed, args.output_dir, args.predictions)

scripts/run_all.slurm

@@ -0,0 +1,30 @@
+#!/bin/bash
+#SBATCH --job-name=qcpinn_cavity_public
+#SBATCH --cpus-per-task=4
+#SBATCH --mem=16G
+#SBATCH --time=16:00:00
+set -euo pipefail
+cd "${SLURM_SUBMIT_DIR:?Submit from the model package root}"
+export OMP_NUM_THREADS=1
+export OPENBLAS_NUM_THREADS=1
+export MKL_NUM_THREADS=1
+export MPLBACKEND=Agg
+
+pids=()
+for seed in {0..9}; do
+  (
+    python -u scripts/train.py --config config/config.yaml --seed "$seed"
+    python -u scripts/inference.py --config config/config.yaml --seed "$seed" --kind best
+    python -u scripts/inference.py --config config/config.yaml --seed "$seed" --kind last --metrics-only
+    python -u scripts/result.py --config config/config.yaml --seed "$seed"
+  ) &
+  pids+=("$!")
+  if ((${#pids[@]} == 4)); then
+    for pid in "${pids[@]}"; do wait "$pid"; done
+    pids=()
+  fi
+done
+for pid in "${pids[@]}"; do wait "$pid"; done
+python -u scripts/result.py --config config/config.yaml --aggregate
+python -u scripts/inference.py --config config/config.yaml
+python -u scripts/result.py --config config/config.yaml

scripts/train.py

@@ -0,0 +1,153 @@
+"""Formal Table 3 training: 20,000 Adam updates, batch64 per constraint."""
+import argparse
+import csv
+import fcntl
+import json
+import random
+import resource
+import sys
+import time
+from pathlib import Path
+
+sys.path.insert(0, str(Path(__file__).resolve().parents[1]))
+import numpy as np
+import torch
+from common import (ROOT, atomic_json, atomic_save, build_model, evaluate, load_checkpoint,
+                    load_config, load_validation, loss_terms, make_sobol_pools,
+                    require_compute_node, resolved, sample_batches, seed_everything,
+                    selection_score, sha256)
+
+
+def train(config, seed, resume=None):
+    allocation = require_compute_node()
+    if seed not in config["training"]["seeds"]:
+        raise ValueError("Seed is not in the formal experiment list")
+    tcfg = config["training"]
+    seed_everything(seed, tcfg["num_threads"])
+    out = resolved(config, "results") / f"seed_{seed}"
+    weights = resolved(config, "weights") / f"seed_{seed}"
+    out.mkdir(parents=True, exist_ok=True)
+    weights.mkdir(parents=True, exist_ok=True)
+    if (out/"history.csv").exists() and resume is None:
+        raise FileExistsError("Existing training history; use --resume to preserve run provenance")
+    model = build_model(config)
+    pools = make_sobol_pools(config, seed)
+    rng = np.random.default_rng(seed+100000)
+    xval, yval, valids = load_validation(config)
+    np.save(out/"validation_indices.npy", valids)
+    optimizer = torch.optim.Adam(model.parameters(), lr=tcfg["learning_rate"],
+                                 betas=tuple(tcfg["betas"]), eps=tcfg["epsilon"],
+                                 weight_decay=tcfg["weight_decay"])
+    sched = tcfg["scheduler"]
+    scheduler = torch.optim.lr_scheduler.ReduceLROnPlateau(optimizer, mode="min", factor=sched["factor"],
+                                                         patience=sched["patience"], threshold=sched["threshold"])
+    start_epoch, best_score, best_epoch = 1, float("inf"), None
+    accumulated_seconds = 0.0
+    if resume:
+        previous = load_checkpoint(resume, model)
+        if previous["config"] != config or previous["seed"] != seed:
+            raise ValueError("Resume config/seed differs from checkpoint")
+        optimizer.load_state_dict(previous["optimizer_state_dict"])
+        scheduler.load_state_dict(previous["scheduler_state_dict"])
+        start_epoch = previous["epoch"]+1
+        best_score, best_epoch = previous["best_score"], previous["best_epoch"]
+        torch.set_rng_state(previous["torch_rng"])
+        random.setstate(previous["python_rng"])
+        rng.bit_generator.state = previous["batch_rng"]
+        accumulated_seconds = previous["elapsed_seconds"]
+    provenance = {"config": config, "seed": seed, "allocation": allocation,
+                  "reference": str(resolved(config, "reference")),
+                  "reference_sha256": sha256(resolved(config, "reference")),
+                  "parameter_count": sum(p.numel() for p in model.parameters()),
+                  "torch_version": torch.__version__, "epoch_semantics": "one Adam step, 64 points per constraint",
+                  "reference_used_for_gradient_training": False,
+                  "validation_indices": str(out/"validation_indices.npy")}
+    atomic_json(out/"run_manifest.json", provenance)
+    print(f"START seed={seed} allocation={json.dumps(allocation)} parameters=923 epochs=20000 batch=64", flush=True)
+    columns = ["epoch", "total", "physics", "lid", "walls", "initial", "lr", "gradient_norm",
+               "step_seconds", "elapsed_seconds", "val_u_l2", "val_v_l2", "val_p_l2", "val_mean_l2"]
+    mode = "a" if resume else "w"
+    start = time.monotonic()
+    with open(out/"history.csv", mode, newline="", buffering=1) as history:
+        writer = csv.DictWriter(history, fieldnames=columns)
+        if not resume: writer.writeheader()
+        for epoch in range(start_epoch, tcfg["epochs"]+1):
+            tick = time.monotonic()
+            model.train()
+            optimizer.zero_grad(set_to_none=True)
+            terms = loss_terms(model, sample_batches(pools, tcfg["batch_size"], rng), config)
+            total = terms["total"]
+            if not torch.isfinite(total): raise FloatingPointError(f"Nonfinite loss at {epoch}")
+            total.backward()
+            grad_norm = torch.nn.utils.clip_grad_norm_(model.parameters(), tcfg["gradient_clip"], error_if_nonfinite=True)
+            optimizer.step()
+            values = {k: float(v.detach()) for k, v in terms.items()}
+            scheduler.step(values["total"])
+            step_seconds = time.monotonic()-tick
+            elapsed = accumulated_seconds+time.monotonic()-start
+            row = {"epoch": epoch, **values, "lr": optimizer.param_groups[0]["lr"],
+                   "gradient_norm": float(grad_norm), "step_seconds": step_seconds, "elapsed_seconds": elapsed}
+            metrics = None
+            if epoch == 1 or epoch % config["validation"]["every"] == 0 or epoch == tcfg["epochs"]:
+                metrics = evaluate(model, xval, yval, config["validation"]["batch_size"])
+                score = selection_score(metrics)
+                row.update({f"val_{k}_l2": metrics[k]["relative_l2"] for k in ("u", "v", "p")})
+                row["val_mean_l2"] = score
+                if score < best_score:
+                    best_score, best_epoch = score, epoch
+                    is_best = True
+                else: is_best = False
+                state = {"model_state_dict": model.state_dict(), "config": config, "seed": seed,
+                         "epoch": epoch, "validation_metrics": metrics, "selection_score": score,
+                         "best_score": best_score, "best_epoch": best_epoch,
+                         "optimizer_state_dict": optimizer.state_dict(), "scheduler_state_dict": scheduler.state_dict(),
+                         "torch_rng": torch.get_rng_state(), "python_rng": random.getstate(),
+                         "batch_rng": rng.bit_generator.state, "elapsed_seconds": elapsed,
+                         "allocation": allocation, "loss": values}
+                if is_best:
+                    atomic_save(weights/"best_model.pt", state)
+                    # Publish the validation-best checkpoint at the requested root path
+                    # throughout concurrent training, with one writer at a time.
+                    root_weights = resolved(config, "weights")
+                    root_meta = resolved(config, "results")/"best_checkpoint_selection.json"
+                    with open(root_weights/".best_model.lock", "a") as lock:
+                        try:
+                            fcntl.flock(lock, fcntl.LOCK_EX | fcntl.LOCK_NB)
+                        except BlockingIOError:
+                            # Another seed is publishing; per-seed checkpoint remains authoritative.
+                            pass
+                        else:
+                            incumbent = json.loads(root_meta.read_text()) if root_meta.exists() else None
+                            if incumbent is None or score < incumbent["selection_score"]:
+                                atomic_save(root_weights/"best_model.pt", state)
+                                atomic_json(root_meta, {"seed": seed, "epoch": epoch, "selection_score": score,
+                                                       "metrics": metrics, "checkpoint": str(root_weights/"best_model.pt")})
+                            fcntl.flock(lock, fcntl.LOCK_UN)
+                if epoch % tcfg["checkpoint_every"] == 0 or epoch == tcfg["epochs"]:
+                    atomic_save(weights/"last_model.pt", state)
+                atomic_json(out/"progress.json", {"status": "training", "epoch": epoch,
+                            "epochs": tcfg["epochs"], "loss": values, "validation_metrics": metrics,
+                            "best_score": best_score, "best_epoch": best_epoch, "allocation": allocation})
+            writer.writerow(row)
+            line = f"seed={seed} epoch={epoch}/20000 loss={values['total']:.8e} phy={values['physics']:.6e} lid={values['lid']:.6e} walls={values['walls']:.6e} ic={values['initial']:.6e} lr={row['lr']:.6e}"
+            if metrics is not None:
+                line += " " + " ".join(f"{k}_relL2={metrics[k]['relative_l2']:.8e} ({metrics[k]['relative_l2_percent']:.5f}%)" for k in ("u", "v", "p"))
+            print(line, flush=True)
+    elapsed = accumulated_seconds + time.monotonic()-start
+    atomic_json(out/"training_summary.json", {"status": "completed", "completed_epochs": tcfg["epochs"],
+                "seed": seed, "best_epoch": best_epoch, "best_validation_mean_l2": best_score,
+                "last_loss": values, "elapsed_seconds": elapsed,
+                "peak_rss_mb": resource.getrusage(resource.RUSAGE_SELF).ru_maxrss/1024,
+                "best_checkpoint": str(weights/"best_model.pt"), "last_checkpoint": str(weights/"last_model.pt"),
+                "allocation": allocation})
+    atomic_json(out/"progress.json", {"status": "completed", "epoch": tcfg["epochs"], "best_epoch": best_epoch})
+    print(f"TRAINING_COMPLETED seed={seed} epochs=20000 best_epoch={best_epoch}", flush=True)
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("--config", default=str(ROOT/"config/config.yaml"))
+    parser.add_argument("--seed", type=int, required=True)
+    parser.add_argument("--resume", type=Path)
+    args = parser.parse_args()
+    train(load_config(args.config), args.seed, args.resume)

weight/best_model.pt

@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:498bcc0dbfc053ca4824b1b6cd2914e85536d0aed1e07ad50d5ef3b4a0b11b87
+size 34274