Repository navigation
fix(solver): apply ns-grid forcing in every PhiFlow RK4 stage - #310
Closed
dionhaefner wants to merge 13 commits into
Closed
dionhaefner wants to merge 13 commits into
dionhaefner wants to merge 13 commits into
Conversation
## Summary `run_timed_trials` makes one warm-up call before each sweep point's timed trials and throws its time away. That call carries what a practitioner pays once: compilation, setup caches, first-touch allocation. This keeps it as `first_call_s` in every timing record (success, and failure after the warm-up), outside `mean`, and adds a "Forward first call" column to the results-page cost table, at the same N as the forward time. Ranking is unchanged. ## Type of change - [ ] New solver backend - [ ] Solver tuning / improvement - [ ] New benchmark domain - [x] Harness / infrastructure change - [ ] Documentation - [ ] Bug fix ## Benchmark label `benchmark:none`. No solver output or existing metric changes; the new field fills in on the next run of each solver. Until then the column shows "—". ## Notes Checked end to end with `mosaic run -p structural-mesh -s FEniCS -e cost/spatial_cost`: | N | `mean` | `first_call_s` | | --- | --- | --- | | 4 | 0.031 s | 18.4 s | | 6 | 0.038 s | 0.10 s | | 32 | 0.42 s | 1.44 s | | 128 | 6.15 s | 13.7 s | Sweep points share a container, so one-off work for all sizes (here the form compilation) lands only on the first point; later points carry only per-size setup. The cost module docstring says so. The table column takes the largest N, so it shows the per-size part. Unit tests cover the warm-up time staying out of the trials and the two failure paths.
## Summary Every pull request opened with a benchmark label shows a failed "Benchmark checks passed" check next to the passing one. Opening with a label fires both `opened` and `labeled`, the concurrency group cancels the first run, and that run's merge gate fails on the cancelled plan. The newer run decides the merge, so the red check is noise, but it reads as a failure to anyone scanning the checks. The gate now passes when its plan was cancelled and a newer run of the workflow exists for the same commit. A run cancelled by hand, with no newer run, still fails, so the required check cannot be satisfied by cancelling. ## Type of change - [ ] New solver backend - [ ] Solver tuning / improvement - [ ] New benchmark domain - [x] Harness / infrastructure change - [ ] Documentation - [ ] Bug fix ## Notes The lookup is one `gh api` call on the workflow's runs filtered by head SHA. Against #288's two runs it finds the newer, successful run for the cancelled one. This pull request only touches `.github/`, so it needs no benchmark label.
…or version (#289) #### Relevant issue or PR None. #### Description of changes PICT's image compiles PISOtorch's CUDA kernels once per entry in `TORCH_CUDA_ARCH_LIST`, and that compile dominates the build. The list was `7.0;7.5;8.0;8.6;9.0`; it is now `7.5;8.0;9.0+PTX`. - 7.5 and 8.0 are native binaries for the two GPUs the benchmark runs on: the T4 in CI and the A100. - A cubin for sm_X.Y also runs on any sm_X.Z with Z ≥ Y, so 8.0 covers 8.6 and 8.9. - `+PTX` embeds compute_90 PTX, which the driver can JIT-compile on GPUs newer than sm_90. Whether the image actually runs on those also depends on the pinned torch build, which I haven't checked. - GPUs older than sm_75 (V100) are no longer supported. The benchmark no longer runs on them. #### Testing done On an A100 box, with the layers before the compile cached in both builds, for an earlier version of this branch with `7.0;8.0;9.0+PTX` (also three architectures): | | main | three architectures | |---|---|---| | image build | 939 s | 610 s | | `PISOtorch*.so` | 56.7 MB | 47.5 MB | - The generated `build.ninja` listed one `sm_` target per architecture plus `compute_90`. - `tests/test_pict_cylinder.py` passed against that image (4 passed). The final list swaps 7.0 for 7.5 and has not been built locally. The `benchmark:solver` run builds it and runs PICT on the T4.
…287) ## Summary PICT's float32 default tolerance for its advection and pressure linear solves is 1e-5. On the periodic box, the residual error in the forward and adjoint solves caps the VJP's agreement with finite differences at about 1e-3, whatever the step size. This passes 1e-6 for both solves in the periodic runs; the cylinder channel is unchanged. ## Type of change - [ ] New solver backend - [ ] Solver tuning / improvement - [ ] New benchmark domain - [ ] Harness / infrastructure change - [ ] Documentation - [x] Bug fix ### For solver tuning - [x] Before/after results compared (below) - [x] No regressions to other solvers ## Benchmark label `benchmark:solver`. ## Notes `mosaic run -s PICT -p ns-grid` (forward, `fd_check`, `analytic_check`, `param_sweep`) on an A100, `main` vs this branch. Every status stays `ok`; forward errors and `analytic_check` agree to three digits, because the spatial discretization dominates there. | ns-grid `fd_check`, median rel. error | `main` | this PR | | --- | --- | --- | | ε = 5 | 8.8e-4 | 8.4e-4 | | ε = 1 | 1.1e-3 | 5.7e-5 | | ε = 0.1 | 1.4e-3 | 9.1e-5 | | ε = 0.01 | 1.1e-3 | 8.8e-4 | The reported error (best ε) goes from 8.8e-4 to 5.7e-5. `param_sweep` mean rel. error at ε = 5 goes from 0.014–0.022 to 0.0013–0.0086 across ν. Without this, PICT's `fd_check` becomes an anomaly once #282 lowers the peer median. Not verified locally: ns-3d-grid, the remaining gradient experiments, and cost. The GPU was shared with other jobs, so timings are not comparable. A tighter tolerance means more solver iterations per step, so expect forward and VJP cost to rise; the `benchmark:solver` run will show by how much.
…293) ## Summary The cost checks (`max_peer_k`) compared each solver's median wall time with the median of all solvers. A serial CPU solver was therefore flagged for being slower than GPU solvers, which says something about the hardware, not the solver, and OpenFOAM carried an explained anomaly for exactly that. Each solver's peers are now the solvers on its device class (`uses_gpu`), the same comparison the results pages make. A class with a single solver gets no peer comparison, as before with fewer than two solvers. ## Effect on the current baseline `mosaic status` on the latest published baseline, `main` against this branch. Only these cells change: | Cell | `main` | This branch | |---|---|---| | ns-grid `cost/spatial_cost` OpenFOAM | anomaly, explained (27× all-solver median) | ok | | ns-grid `cost/temporal_cost` OpenFOAM | anomaly, explained (20×) | ok | | ns-grid `cost/spatial_cost` PhiFlow | ok | anomaly (42× GPU median: 4.0 s against 0.10 s) | | ns-grid `cost/temporal_cost` PhiFlow | ok | anomaly (21×) | | ns-grid `cost/vjp_cost/by_N` PhiFlow | ok | anomaly (21×) | The PhiFlow anomalies are real: against the other GPU solvers, PhiFlow's 2D runs are 20–40× slower. The CPU solvers in the pooled median hid this. Earlier measurements found PhiFlow's 2D case slower on the GPU than on the CPU, so the next step is to find out why, not to explain it away. ns-grid's campaign score goes from 0.962 to 0.958. ## Type of change - [ ] New solver backend - [ ] Solver tuning / improvement - [ ] New benchmark domain - [x] Harness / infrastructure change - [ ] Documentation - [ ] Bug fix ## Benchmark label `benchmark:none`: no solver output changes, only how the status checks classify existing results. ## Notes - Removes `OPENFOAM_SERIAL_IMPLICIT_COST`, which no longer applies. - Four new tests in `tests/core/test_status.py`. - Merging #288 (PhiFlow with RK4 in the periodic box, 2.6–3× the cost) will make the PhiFlow anomalies larger.
…d XLB (#291) ## Summary jax-cfd, PhiFlow and XLB ran with `--xla_gpu_autotune_level=0 --xla_gpu_enable_cublaslt=false`, which switches off two optimizations XLA enables by default. The reasons on record are all from V100 GPUs: an autotuner `DEVICE_TYPE_INVALID` error (XLB), an out-of-memory error in PhiFlow's 3D recovery on 16 GiB cards, and a cuBLASLt initialization crash on a busy GPU. The benchmark now runs on T4s in CI and on A100s, so this restores the framework defaults, following the tuning checklist: a default optimization is only switched off for a reason that reproduces on the benchmark hardware. ## Type of change - [ ] New solver backend - [x] Solver tuning / improvement - [ ] New benchmark domain - [ ] Harness / infrastructure change - [ ] Documentation - [ ] Bug fix ### For solver tuning - [ ] Before/after `mosaic status --format json` snapshots compared - [ ] No regressions to other solvers ## Benchmark label `benchmark:solver`. The run is the test: it shows whether any of the three solvers fails on the T4, and how cost changes. Autotuning compiles extra kernels in the untimed warm-up call, so steady-state cost should not get worse. ## Notes - If a failure reproduces on the T4, the flag goes back for that solver, with the T4 as the recorded reason. PhiFlow is the most likely case, since the T4 also has 16 GiB. - Not yet checked on an A100.
…ox and the channel (#288) ## Summary The PhiFlow tesseract advects velocity with central differences and forward Euler, a pairing PhiFlow's documentation never uses. This switches to the documented schemes: central differences under `fluid.incompressible_rk4` in the periodic boxes ("Higher-Order Simulations" in PhiFlow's fluid guide), and the default semi-Lagrangian advection with operator splitting in the cylinder channel, since `incompressible_rk4` takes no obstacles. **Periodic forward and VJP cost rise ×2.4–3.9 and ×3.1–4.4**: RK4 makes four pressure projections per step instead of one. ## Type of change - [ ] New solver backend - [x] Solver tuning / improvement - [ ] New benchmark domain - [ ] Harness / infrastructure change - [ ] Documentation - [ ] Bug fix ### For solver tuning - [x] Before/after results compared (below) - [x] No regressions to other solvers ## Benchmark label `benchmark:solver`. ## Notes `mosaic run -s PhiFlow` on an A100, `main` vs this branch; cost timed with one job on the GPU. | this PR / `main` | ns-grid | ns-3d-grid | | --- | --- | --- | | forward | ×2.6–3.0 (N = 256: 14.5 → 41.0 s) | ×2.4–3.9 (N = 64: 3.96 → 14.1 s) | | VJP | ×3.1–3.4 (N = 256: 29.5 → 96.9 s) | ×3.3–4.4 (N = 64: 8.1 → 35.9 s) | | | `main` | this PR | | --- | --- | --- | | ns-grid `baseline` error | 5.94e-3 | 5.75e-3 | | ns-grid `cylinder` error (Cd vs DFG) | 0.0538 | 0.0111 | | ns-grid `fd_check` rel. error | 9.1e-6 | 1.3e-5 | | ns-3d-grid `baseline` error | 0.041 | 0.0398 | | ns-3d-grid `fd_check` rel. error | 3.2e-5 | 1.2e-5 | | ns-3d-grid `horizon_sweep_limits` | NaN at 160 steps | NaN at 160 steps | Periodic accuracy barely moves; the spatial error dominates. The reasons are attribution and stability: forward Euler with central advection is stable only for roughly dt ≤ 2ν/|u|², which vanishes as ν → 0. On a 2-D Taylor–Green vortex (N = 32, ν = 0, dt = 0.05), `main` goes NaN at 800 steps; this branch decays monotonically to 2000. Today's checks stay inside that limit or run too few steps to show it. The 3-D horizon failure is unchanged: both schemes gain energy and fail at the same step, at a CFL number of 0.16. Not rerun here: `drag_opt` (semi-Lagrangian everywhere gave 47.6% in 65 iterations vs 46.8% in 57). `incompressible_rk4` in PhiFlow 3.4.0 is third order (stage 4 uses the second slope); that does not show here.
… caches (#292) ## Summary The setup caches in JAX-FEM (structural and thermal), TopOpt.jl and INS.jl kept one entry per mesh a sweep visited, so a container's memory grew with every sweep point and the largest problem size ran with less headroom than a single setup needs. They now hold the current setup only, as the FEniCS and Firedrake caches already do. This is item 6 of the harness contract for runtime adaptations: "a cache holds the current setup only". ## Type of change - [ ] New solver backend - [x] Solver tuning / improvement - [ ] New benchmark domain - [ ] Harness / infrastructure change - [ ] Documentation - [ ] Bug fix ## Benchmark label `benchmark:solver`. Outputs should be unchanged, because the timed calls at each sweep point reuse the setup built in the warm-up call. The recorded memory metrics cannot show the effect: `ram_peak_mib` is the harness process's RSS, not the container's, and `vram_peak_mib` is a whole-GPU delta per point. ## Notes An earlier implementation of the same four changes, made on the A100 sandbox, left every output bit-identical over the forward, gradient and cost suites. This branch is a fresh implementation; the benchmark run checks it.
… once per mesh (#286) ## Summary Both JAX-FEM tesseracts solve with `umfpack_solver` (SciPy `spsolve` on the host), copied from the topology-optimization demo. JAX-FEM's default `jax_solver` (Jacobi-preconditioned BiCGSTAB on the GPU) is 5 to 8 times faster on large meshes, but it costs about 0.5 s per forward and 1 s per VJP even on tiny ones. That overhead is recompilation: `jax_solve` runs `bicgstab` eagerly, so it compiles 2 XLA programs per forward and 4 per VJP. This passes the same solve, jitted once per matrix size, through JAX-FEM's `custom_solver` option. After the first call, compilations drop to 0. ## Type of change - [ ] New solver backend - [x] Solver tuning / improvement - [ ] New benchmark domain - [ ] Harness / infrastructure change - [ ] Documentation - [ ] Bug fix ### For solver tuning - [x] Before/after results compared (below) - [x] No regressions to other solvers ## Benchmark label `benchmark:solver`. ## Notes Draft until the timing below is repeated on a quiet GPU. Accuracy (all suites, A100): every metric is identical to umfpack to the recorded digits, except the thermal `conductivity_recovery_bfgs` final ratio (2.17e-8 → 2.4e-8, iterative solve at tol 1e-10). Cost, mean of 3 trials in seconds. These were measured on a GPU shared with other runs: | | N | umfpack (`main`) | `jax_solver` | this PR | | --- | --- | --- | --- | --- | | structural `spatial_cost` | 4 / 128 | 0.26 / 14.0 | 0.56 / 2.56 | 0.30 / 2.46 | | structural `vjp_cost/by_N` | 4 / 128 | 0.48 / 24.5 | 1.13 / 5.62 | 0.62 / 4.50 | | thermal `spatial_cost` | 16 / 1024 | 0.23 / 69.6 | 0.57 / 7.0 | 0.29 / 8.4 | | thermal `vjp_cost/by_N` | 16 / 1024 | 0.58 / 137 | 1.2 / 14.3 | 0.64 / 14.0 | At small N a call is still 15 to 30% slower than umfpack. Wall time of the small-mesh optimizations: thermal recovery 1026 → 1012 s, structural `topopt` 2591 → 4241 s. The two runs overlapped, so I don't trust the topopt number yet.
…e rule (#270) ## Summary TopOpt.jl's VJP read the undocumented `Compliance.grad` field after a forward call. It now takes the gradient through Zygote, as TopOpt.jl documents (`Zygote.gradient(r -> comp(PseudoDensities(r)), rho)`). TopOpt.jl 0.11.0's `rrule` for `Compliance` returns `copy(comp.grad)` scaled by the cotangent, so the gradient is the same vector. `topopt_forward` now returns only the compliance; the gradient and displacement it also returned were unused apart from the VJP. Stacked on #260 (the file moves there); retarget to `main` once #260 is merged. ## Type of change - [ ] New solver backend - [ ] Solver tuning / improvement - [ ] New benchmark domain - [ ] Harness / infrastructure change - [ ] Documentation - [ ] Bug fix Cleanup, no result change intended. ## Benchmark label `benchmark:solver`, so CI confirms TopOpt.jl results are unchanged. ## Notes Verified locally (CPU): - In the image, `comp.grad` and the Zygote gradient are bit-identical on cantilevers at nx = 8 and 16 with random densities. - `gradient/fd_check`, `gradient/param_sweep` and `optimization/topopt` give identical metrics with images from #260 and from this branch (final compliance 0.0010033302241936326 after 2500 iterations in both). - `cost/vjp_cost` (`by_N`, `by_steps`), three runs per image: VJP time is the same within run-to-run spread, e.g. 0.66–0.72 s with `comp.grad` and 0.65–0.99 s with Zygote at nx = 128 (the 0.99 s run was an outlier; the other two gave 0.65 and 0.67 s).
#### Relevant issue or PR Found while adding the Kolmogorov solver-in-the-loop experiment (#274), where PICT fails the finite-difference check of the training gradient. #### Description of changes **Deviation from PICT.** PICT's linear solves stop once the RMS of their residual falls below an absolute tolerance, 1e-5 in float32 by default. The adjoint solves of its backward pass use the same tolerance. When the cotangent is small, the adjoint residual starts close to or below that tolerance, so the adjoint solves stop early. PICT's VJP is then not linear in the cotangent, and the gradient depends on how the loss is scaled. Single calls with O(1) cotangents, which is what the gradient suite checks, are unaffected. The problem shows up when a scalar loss sits at the end of a chain of calls. On the Kolmogorov solver-in-the-loop loss (16 chained calls, cotangents of about 1e-4 per entry), the gradient with respect to the state after interval k disagrees with finite differences by: | k (calls left to the loss) | 15 (1) | 13 (3) | 8 (8) | 0 (16) | |---|---|---|---|---| | relative error, as is | 0.04% | 3.7% | 10% | 73% | | loss × 1e4 | — | 0.014% | 0.023% | 0.015% | The tesseract's VJP now scales the cotangents to a largest entry of 1 before PICT's backward pass and scales the gradients back. This is exact for a linear VJP and makes the adjoint tolerance relative to the cotangent. The code comment and the tesseract's description (rendered on the docs solvers page) both state the deviation. An alternative is a relative convergence criterion in PISOtorch. PICT's criteria are all absolute (`ABS_MAX`, `ABS_MEAN`, `ABS_SUM`, `NORM2`, `NORM2_NORMALIZED`), so that would mean patching PICT itself. #### Testing done - New `tests/test_pict_vjp.py`: the VJP for a cotangent scaled by 1e-5, divided by 1e-5, matches the VJP for the unscaled cotangent to within 1e-3. The test fails on main's image (4.3%) and passes with this change. `tests/test_pict_cylinder.py` passes too. - A/B of PICT's gradient and optimization suites on ns-grid and ns-3d-grid, main's image against this branch's (local A100 run, same harness): | | main | this PR | |---|---|---| | ns-grid `fd_check`, rel. error at ε = 1 / 0.1 / 0.01 | 1.1e-3 / 1.4e-3 / 1.1e-3 | 2.7e-4 / 3.2e-4 / 7.2e-4 | | ns-3d-grid `fd_check`, same | 1.1e-4 / 5.0e-4 / 4.6e-3 | 1.4e-4 / 3.7e-4 / 4.5e-3 | | ns-grid `analytic_check` error, ν = 1e-3 … 5e-2 | 0.0043 … 0.0198 | 0.0044 … 0.0198 | | ns-3d-grid recovery, final IC error (3 trials) | 0.0013, 0.0015, 0.0021 | 0.0022, 0.0028, 0.0040 | | ns-3d-grid recovery, final loss | 3.5e-10 … 6.4e-10 | 5.7e-10 … 8.5e-10 | In the other gradient experiments (horizon sweeps, Jacobian spectra), gradient norms and condition numbers change by 6% at most, and every status stays the same. PICT is excluded from `drag_opt`, so recovery is the only optimization it runs. Recovery gets slightly worse: both images drive the loss to the same float32 floor, but the IC error ends about twice as large with this change. A repeat of both arms reproduced every trial to two digits, so it is not run-to-run noise. Both are far below the failure threshold of 2.0 (from an initial error of 1.0). I don't have an explanation that I've verified; the loss is very flat there (the Jacobian's condition number is about 7e7), so where L-BFGS stops in the directions the loss barely sees could depend on details of the gradient.
## Summary Both deal.II tesseracts launched the C++ solver once per call, so every call paid about 0.2 s for process start-up and for loading deal.II. On the small meshes that is almost the whole cost. The binaries gain a `--serve` mode that reads one `input.json` path per line from stdin and answers `ok` or `error: <message>`. Each tesseract keeps one such process alive and restarts it if it exits. Each request still builds its own mesh and linear system, so nothing that depends on an input carries over between calls. This is the persistent worker that the harness contract counts as a setup cache (policy 2, work item 3). Single-shot use (`heat_solver input.json`) is unchanged. ## Type of change - [ ] New solver backend - [x] Solver tuning / improvement - [ ] New benchmark domain - [ ] Harness / infrastructure change - [ ] Documentation - [ ] Bug fix ### For solver tuning - [x] Before/after results compared (below) - [x] No regressions to other solvers ## Benchmark label `benchmark:solver`. ## Notes Local A/B, every suite on both problems. Non-cost metrics are identical in all 47 cells. Cost, mean of 3 trials: | | N | `main` | this PR | | --- | --- | --- | --- | | structural `spatial_cost` | 4 / 32 / 128 | 0.206 / 0.332 / 2.90 s | 0.009 / 0.131 / 2.66 s | | structural `temporal_cost` | | 0.217 s | 0.016 s | | thermal `spatial_cost` | 16 / 128 / 1024 | 0.220 / 0.650 / 52.9 s | 0.013 / 0.411 / 57.3 s | | thermal `temporal_cost` | | 0.300 s | 0.101 s | Image builds shared the CPU during part of both runs, so the large-mesh ratios (−27% to +8%) are indicative only. The small-mesh gain is far larger than that noise. This PR's benchmark run is the clean measurement.
… the RK4 momentum main's PhiFlow (#288) integrates the periodic box with incompressible_rk4, so the body force and the linear drag move into its momentum rate and every RK4 stage sees them. PICT keeps both sides of its conflict.
Contributor
Author
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to subscribe to this conversation on GitHub.
Already have an account?
Sign in.
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Merges
mainintofeat/ns-grid-forcing(#267) and resolves the three conflicts. The one that changes numerics is PhiFlow. Since #288,mainintegrates the periodic box withincompressible_rk4, so the body force and the linear drag move into its momentum rate. Every RK4 stage now sees them, and the drag acts on that stage's velocity.Type of change
What changes
momentum()underincompressible_rk4. The docstring now says "added to the momentum rate".prep_fnvelocity source andmain's tolerance comment.main, which removed the KBC collision selector. The forcing comment no longer mentions KBC.Verification
tests/test_tmp_forcing.pypasses on images built from this branch: the laminar Kolmogorov transient, zero force and drag change nothing, and the forced VJP matches central differences, for PhiFlow and PICT (6 passed). The OpenFOAM test is skipped because its image isn't built here. The XLB forcing path has no test in that file and was not run.Benchmark label
benchmark:solver: PhiFlow's forced time stepping changes.