diff --git a/CMakeLists.txt b/CMakeLists.txt index b652aa6..5fe997d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,7 +2,14 @@ cmake_minimum_required(VERSION 3.14) if(POLICY CMP0135) cmake_policy(SET CMP0135 NEW) endif() -project(my-fantastic-simulation-code LANGUAGES C CXX VERSION 0.1) +project(qew LANGUAGES C CXX VERSION 0.1) + +# Bandwidth-bound HPC code: an accidental no-flags build (-O0) is never what +# anyone wants, so default to Release when the user does not choose. +if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) + set(CMAKE_BUILD_TYPE Release CACHE STRING "Build type" FORCE) + message(STATUS "CMAKE_BUILD_TYPE not set, defaulting to Release") +endif() # Kokkos 5.x requires C++20 or newer. set(CMAKE_CXX_STANDARD 20) diff --git a/docs/fix_plan.md b/docs/fix_plan.md new file mode 100644 index 0000000..e74d510 --- /dev/null +++ b/docs/fix_plan.md @@ -0,0 +1,176 @@ +# Implementation plan for `optimization_notes.md` + +> **Status (June 2026): all phases 0–4 are implemented on this branch.** +> Phase 1 was settled by the paper (arXiv:2410.21838 eq. (10): +> `lambda_p/xi = (Gamma/U0)^{2/3}`, varied by scanning Gamma at fixed +> disorder): all drivers now take the dimensionless ratio as +> `pinning_length_over_xi` and share `Gamma = A (lambda_p/xi)^{3/2}`, with +> per-driver mapping tests; line tensions are numerically unchanged. +> Measured on the benchmark (Serial): Phase 2 took the end-to-end +> preconditioned L-BFGS relaxation (n = 4096, arclength) from 5.73 s / 12893 +> objective evals to 4.77 s / 10424 at an identical 2470 iterations; Phase 3 +> took it to 3.55 s (38% total) and the hot kernels (n = 262144) from +> grad 33.5 / dyn 52.4 ms to 26.5 / 34.4 ms, with the L-BFGS trajectory +> bit-identical throughout. Phase 4's RK change removes one of two +> reductions per sweep; the Langevin RNG chunking is serial-neutral and its +> GPU benefit still needs the CUDA runner (see below). Phase 5 remains +> won't-fix unless profiling says otherwise. + +Companion to [`optimization_notes.md`](optimization_notes.md); item numbers +below refer to that document. Work is grouped into five phases that are +independently land-able, ordered so that correctness fixes go first and each +perf change is benchmarked in isolation. Every phase ends with the full +`ctest` suite on a Serial **and** an OpenMP build (the threaded run doubles as +the race check), plus a before/after run of `executables/benchmark` for the +perf phases. + +Ground rules: + +- The gradient-consistency test and the golden-file tests are the correctness + gates; they must pass unchanged (same tolerances) in every phase. +- Tests that assert `n_func_evals` counts will be updated *deliberately* in + Phase 2, in the same commit as the behaviour change, with the new expected + count derived by hand. +- One commit (or small commit series) per phase, so a regression bisects to a + phase. + +## Phase 0 — Safety & build hygiene (items 1, 3, 4) + +Small, riskless, do first. + +1. **Odd-`n` guard** (`rosso_krauth.h`): in `rk_block()`, before launching + kernels, `throw std::invalid_argument` (host code) when `n % 2 != 0`, + citing the red-black neighbour conflict. A 3-colour sweep for odd `n` is + *not* worth the complexity until someone needs odd lines. Add a death/throw + test. +2. **Integer ramp counter** (`qew_runs.h::depinning_ramp`): iterate + `for (int k = 0; k * f_step <= f_max; ++k)` with `f = k * f_step`. The + default `f_step = 1.0` produces bit-identical forces, so no test or HDF5 + key changes. +3. **CMake**: default `CMAKE_BUILD_TYPE=Release` when unset (with a STATUS + message); rename the project to `qew`. Leave `hello.cpp` / + `kokkos_smoketest.cpp` in place for now — they are referenced by tests and + removing them is cosmetic churn; fold into a later cleanup if desired. + +## Phase 1 — Decision point: line-tension convention (item 2) + +Blocked on information outside the repo: the `../Overleaf/qEW/qew_*.py` +sources. **Needs Lars to confirm** which convention each Python driver uses. +Then, in one commit: + +- Rename the field in whichever config struct holds the ratio to + `pinning_length_over_xi` (or convert at the call site so all three drivers + use one formula). +- Add a tiny unit test per driver pinning the `pinning_length → + line_tension` mapping, so the convention can never silently diverge again. + +Nothing later depends on this; it can land any time. + +## Phase 2 — L-BFGS overhead (items 5–7) + +All inside `lbfgs.h`, in three commits, benchmarking `qew_static`-shaped runs +(preconditioned, n = 4096) after each: + +1. **Carry the objective forward** (item 5). Hoist one + `phi = objective(h, ...)` above the iteration loop; on line-search + acceptance set `phi = phi_accepted`; delete the per-iteration `phi0` + recompute. The line-search-failure `continue` path restores `h = x0`, so + the carried value stays valid there too. Arithmetic is unchanged → + golden L-BFGS test passes bit-for-bit; only `n_func_evals` assertions in + `test_lbfgs.cpp` need recounting. +2. **Fuse trial step into the objective** (item 6). Add a local helper + `trial_objective(x0, d, alpha, p, noise)` — the `objective()` reduce with + `h(i)` replaced by `x0(i) + alpha * d(i)` (and the same for the `ip` + neighbour, which is the identical FP expression, so results stay + bit-identical). Commit `h` with a single `parallel_for` only on + acceptance; the `deep_copy(x0, h)` save and `deep_copy(h, x0)` undo both + disappear. Keep `objective()` itself untouched (public API, used by + tests/benchmark). +3. **Reduction fusion + history swap** (item 7). + - One `parallel_reduce` with three reducers computing `sy`, `ss`, `yy` + while also writing `s_tmp`/`y_tmp` (replaces one `parallel_for` + three + reduces). + - `std::swap(S[row], s_tmp)` / `std::swap(Y[row], y_tmp)` instead of + `deep_copy` (safe: the FFT plans are bound to `q`/`d` only). + - Fold `q -= a*y` into the subsequent dot where profiling on a GPU build + shows the sync count matters; otherwise skip — this is the only + speculative sub-item, and it slightly obfuscates the two-loop recursion. + +## Phase 3 — Hot-kernel work reduction (items 8, 9, 10) + +The order matters: land 10 (mechanical, bit-identical) first so the bigger +rewrites diff cleanly against it. + +1. **Modulo removal** (item 10, `filtered_noise.h::sample`, + `qew_model.h`, `rosso_krauth.h`): replace `(i±1) % n` stencil wraps with + compare-and-adjust; in `sample()`, clamp the `gx == nx` rounding edge case + explicitly (`if (ix >= nx) { ix -= nx; gx -= nx; }` before computing + `fx`), then build the four neighbour indices with conditional wraps. + Same arithmetic → golden noise-interpolation test must pass bit-for-bit. +2. **Value-only / dy-only sampling** (item 9). Add to `DeviceNoise`: + `sample_value(x, y)` (no `dw`, no derivative accumulators) and keep + `sample()` as-is. In `qew_model.h`, route calls through C++20 + `requires`-detected helpers (`noise_value(noise, x, y)` falls back to + `noise.sample(x, y).v` for the analytic test stubs, so no test functor + has to change). Split `bspline_weights` into weights-only and + weights+derivatives variants. `objective()` switches to the value path. +3. **Edge-pass gradient** (item 8). Introduce + `struct GradientWorkspace { View1D edge; }` (sized lazily to `n`). + `gradient(h, p, noise, grad, ws)` becomes two kernels: + - edge kernel: for edge `(i, i+1)` store + `edge(i) = 0.5 * dv_dy(midpoint)` — plus, for the arclength model, the + normalized slope `dh/sqrt(1+dh²)`, packed into a second view or a 2-D + workspace, since both adjacent sites reuse it; + - combine kernel: `grad(i) = dx * (line_force(from edge slopes or h) + + edge(im) + edge(i) - f)`. + Keep a back-compat `gradient(h, p, noise, grad)` overload that allocates a + workspace internally (used by tests); `lbfgs_minimize`, `fire_minimize`, + and `LangevinDynamics` own a persistent workspace. Floating-point sums + re-associate slightly (`0.5*(sl+sr)` becomes `0.5*sl + 0.5*sr`), so the + gradient-consistency test (tolerance-based) is the gate, and the golden + L-BFGS comparison tolerance may need a documented epsilon bump — verify + it stays within the existing tolerance first. + `rk_site_velocity` is intentionally untouched (single-site re-evaluation + cannot share edges). + +Acceptance: benchmark `grad` ms/eval roughly halves on OpenMP; `obj` +improves measurably from item 9; no test tolerance loosened beyond what is +justified in the commit message. + +## Phase 4 — Rosso–Krauth & dynamics (items 11, 12) + +1. **Drop the per-sweep mean reduction** (item 11): track + `cum_adv += total_adv` and compute `res.mean_h = mean_h0 + cum_adv / n`. + Bit-identical apart from summation order of the mean; the RK tests + (washboard `f_c`, runaway) are tolerance-based and unaffected. +2. **Langevin RNG** (item 12): measure first. On the OpenMP build the pool + is per-thread and fine; the concern is GPU lock contention. Run the + benchmark's `dyn` timing on the CUDA CI runner; only if `dyn − grad` + overhead is significant, switch to a pre-generated normal-deviate view + filled by a dedicated kernel where each thread owns one pool state for the + whole row-chunk (avoids per-element lock churn), fused with the update + step. The equipartition test gates correctness; exact trajectories may + change (different draw order) — that test is statistical, so no fixture + updates expected. + +## Phase 5 — Optional cleanups (items 13–15) + +FIRE reduction fusion, host-side noise-setup loops, `to_host_vec` memcpy. +Do only if profiling ever shows them; otherwise close as won't-fix and note +that in `optimization_notes.md`. + +## Sequencing summary + +| Phase | Contents | Risk | Blocked on | +|---|---|---|---| +| 0 | odd-`n` guard, ramp counter, CMake | none | — | +| 1 | line-tension convention | none | Python sources (Lars) | +| 2 | L-BFGS: carry φ, fused trial, fused reduces, view swap | low | — | +| 3 | modulo removal → sample variants → edge-pass gradient | medium | — | +| 4 | RK mean tracking; Langevin RNG (measure-first) | low | GPU runner for 12 | +| 5 | FIRE / setup / copy cleanups | — | profiling evidence | + +Phases 0, 2, 3, 4 are mutually independent except that Phase 3's edge-pass +should rebase on Phase 2's line-search changes (both touch the `lbfgs.h` +call sites). Suggested landing order: 0 → 2 → 3 → 4, with 1 whenever the +Python check happens. diff --git a/docs/optimization_notes.md b/docs/optimization_notes.md new file mode 100644 index 0000000..1136e92 --- /dev/null +++ b/docs/optimization_notes.md @@ -0,0 +1,159 @@ +# Optimization potential & known issues + +Technical notes from a code review of the qEW port (June 2026). Items are +ordered by expected impact. + +> **Status:** items 1–12 are addressed on this branch — see `fix_plan.md` +> for the phase breakdown and measured numbers. Item 12's GPU benefit +> (RNG lock relief) still needs verification on the CUDA runner; the minor +> items 13–15 are won't-fix unless profiling says otherwise. File references point at the current `main`-line +sources. Nothing here changes physics; everything preserves the documented +analytic consistency between `objective()` and `gradient()` and the golden-file +contracts. + +The hot path, per the benchmark (`executables/benchmark`), is the bicubic noise +gather inside `gradient()` / `objective()` — memory-bandwidth-bound (16 +coefficient loads + ~70 FLOPs per `DeviceNoise::sample()`). Most of the wins +below attack either (a) the number of `sample()` calls, or (b) per-iteration +kernel-launch / device-sync overhead in the minimizers. + +## Bugs / correctness hazards + +1. **Red–black sweep races for odd `n`** — `src/rosso_krauth.h:147` + (`if ((i & 1) != color) return;`). On a periodic ring with odd `n`, sites + `0` and `n-1` have the same parity but are nearest neighbours, so the same + colour pass updates both while each reads the other's `h`. All current + configs use even `nx`, so this is latent — but nothing enforces it. + *Fix:* assert `n % 2 == 0` in `rk_block()` (cheap), or fall back to a + 3-colour sweep for odd `n`. + +2. **Inconsistent pinning-length → line-tension map across run drivers** — + `src/qew_runs.h`: `static_sweep()` and `depinning_ramp()` used + `line_tension = pinning_length^1.5 * amplitude`, while `dynamic_anneal()` + used `(pinning_length / xi)^1.5 * amplitude` — the same struct-field name + meant two different physical quantities. **Resolved:** the paper + (arXiv:2410.21838, eq. (10)) controls the model through the dimensionless + Larkin ratio `lambda_p/xi = (Gamma/U0)^{2/3}`, so all drivers now take a + `pinning_length_over_xi` field and share `Gamma = A * (lambda_p/xi)^{3/2}`; + each driver's mapping is pinned by a unit test. HDF5 key strings are + unchanged and their values carry the ratio (already true for static and + depinning output; dynamic output previously recorded the physical length). + +3. **Floating-point loop counter in the force ramp** — + `src/qew_runs.h:214`: `for (real_t f = 0; f <= c.f_max; f += c.f_step)`. + Exact for the default `f_step = 1.0`, but for a non-representable step + (e.g. `0.1`) accumulated rounding can drop or shift the final step and + perturb the HDF5 key strings (`driving_force=%.4g`). *Fix:* integer index, + `f = k * f_step`. + +4. **No default `CMAKE_BUILD_TYPE`** — root `CMakeLists.txt` never sets a + fallback, so configuring without `-DCMAKE_BUILD_TYPE=Release` silently + produces an unoptimized (`-O0`) build of a bandwidth-bound HPC code. *Fix:* + default to `Release` when unset (warn). Cosmetic: the project is still + named `my-fantastic-simulation-code` (line 5), and `hello.cpp` / + `kokkos_smoketest.cpp` scaffolding is still compiled into the production + `lib` target. + +## Optimization potential + +### Minimizer (L-BFGS, `src/lbfgs.h`) + +5. **Redundant objective evaluation every iteration** — line 200: + `phi0 = objective(h, ...)` is recomputed at the start of each line search, + but `h` is exactly the point accepted by the previous iteration's line + search, whose `phi` is already known. Carrying the accepted value forward + removes one full objective (one `sample()` per site + reduction) per + iteration — when the first trial step is accepted (the common, + preconditioned case) this *halves* the objective-evaluation cost. + +6. **Fuse the trial step into the objective reduce** — each Armijo trial does + a `parallel_for` (write `h = x0 + a*d`, line 212) plus a `parallel_reduce` + (objective). Evaluating the objective directly at `x0(i) + a*d(i)` inside + the reduce, committing `h` only on acceptance, saves one kernel launch and + `n` stores per trial and removes the `deep_copy(h, x0)` undo path + (line 224) and the `deep_copy(x0, h)` save (line 199). + +7. **Fewer reductions / device syncs per iteration** — + - `sy`, `ss`, `yy` (lines 245–247) are three separate `parallel_reduce`s + over the same two vectors; fuse into one reduce with a 3-component + reducer. The `lbfgs_sy` `parallel_for` (line 240) can be folded in too. + - The two-loop recursion issues `2m` blocking dot products per iteration, + each a device→host sync; for the production `n = 4096` a GPU run is + latency- rather than bandwidth-limited here. The `axpy` can be fused + with the next loop's dot (`q -= a*y` and `s·q` in one kernel). + - History update: `deep_copy(S[row], s_tmp)` / `deep_copy(Y[row], y_tmp)` + (lines 256–257) copy two full vectors per iteration; `std::swap`-ing the + `View`s into the slot is free. (The FFT plans are bound to `q`/`d`, not + `s_tmp`/`y_tmp`, so swapping is safe.) + +### Force kernel (`src/qew_model.h`) + +8. **Each noise edge is sampled twice in `gradient()`** — site `i` samples the + midpoints of edges `(i-1,i)` and `(i,i+1)` (lines 105–106), so every edge + is evaluated by both adjacent sites: `2n` bicubic gathers for `n` distinct + edges. A two-pass scheme — edge kernel writing `dv_dy` (and the arclength + slope term) per edge into a scratch view, then a cheap combine kernel — + does `n` gathers + `~3n` doubles of scratch traffic instead of `2n` + gathers (16 loads each). Expect close to 2× on the dominant kernel. + (`objective()` already samples each edge exactly once; `rk_site_velocity` + re-evaluates at varying `h` and cannot share.) + +9. **`sample()` computes derivatives the caller discards** — + `objective()` uses only `s.v`, yet `DeviceNoise::sample()` + (`src/filtered_noise.h:67`) always accumulates `dv_dgx` *and* `dv_dgy` and + `bspline_weights()` always fills `dw[4]`. `gradient()` never uses `dv_dx` + anywhere in the codebase. Splitting into value-only / value+dy variants + (template tag or separate methods) trims ~⅓–½ of the FLOPs in the 4×4 + inner loop. Gain is modest where memory-bound, but the Armijo line search + is objective-only and benefits directly. + +10. **Integer modulo in the inner loops** — each `sample()` does 8 `%` + operations for the neighbour wrap (`src/filtered_noise.h:79–80`), plus the + stencils' `(i±1) % n`; integer division is expensive on GPUs. After the + positive wrap, `ix ∈ [0, nx]` (the `gx == nx` rounding edge case is + currently absorbed by `ix % nx`), so conditional add/subtract wraps are + sufficient — keep an explicit clamp for that edge case. Power-of-two + grids could use bit masks. + +### Rosso–Krauth (`src/rosso_krauth.h`) + +11. **Redundant mean-`h` reduction per sweep** — `rk_block()` runs two + reductions per sweep (lines 161–167): the advance sum and a fresh mean of + `h`. Since every move is forward, `mean_h = mean_h0 + Σ_sweeps + total_adv / n`; tracking the cumulative advance eliminates the second + reduction. Sweeps run up to `max_sweeps = 2·10⁵`, so this halves the + per-sweep reduction/sync count. + +### Langevin dynamics (`src/dynamics.h`) + +12. **RNG pool acquisition per site per step** — `pool.get_state()` / + `free_state()` inside the `parallel_for` (lines 48–50) acquires a locked + state per site per step. With the benchmark's `nx = 262144` on a GPU the + default-sized `Random_XorShift64_Pool` serializes on lock contention. + Alternatives: size the pool to the launch width, or pre-fill a normal + deviate view per step with `Kokkos::fill_random` and fuse the update. + Also, the `gradient()` call inside `step()` inherits item 8's 2× edge + duplication. + +### Minor / one-time + +13. **FIRE** (`src/fire.h`) — `nv2`/`nf2` are two separate reductions + (lines 103–108) and the v-step (line 79) could fuse with the power + reduction (line 85). Low priority: FIRE is the fallback minimizer. +14. **Noise setup** (`src/filtered_noise.cpp`) — phase generation, mean and + variance are serial host loops over `nx·ny`; `build_coefficients()` stages + through three complex 2-D buffers + host mirrors. One-time cost + (`setup ≈` seconds in the benchmark), only worth touching if many seeds + are ever swept in one run. +15. **`to_host_vec`** (`src/qew_runs.h:28`) — extra element-wise copy from the + mirror into a `std::vector`; could `memcpy` from the mirror when layouts + match. Negligible at current snapshot rates. + +## Suggested order of attack + +Items 5–7 (L-BFGS overhead) and 8 (edge-duplication) are independent, +low-risk, and cover the production `qew_static` runtime; the +gradient-consistency and golden-file tests pin their correctness. Item 11 is +a five-line change covering `qew_depinning`. Item 1 should be fixed (or +asserted away) before anyone runs an odd-`n` line. Items 2 and 4 are +correctness-of-use issues rather than speed, and are cheap insurance. diff --git a/executables/benchmark/main.cpp b/executables/benchmark/main.cpp index 516e34a..0a36599 100644 --- a/executables/benchmark/main.cpp +++ b/executables/benchmark/main.cpp @@ -10,6 +10,7 @@ #include "dynamics.h" #include "filtered_noise.h" +#include "lbfgs.h" #include "qew_model.h" #include "qew_types.h" @@ -69,6 +70,35 @@ int main(int argc, char *argv[]) { "dyn=%7.3f ms/ea\n", nthreads, t_setup, 1e3 * t_grad / reps, 1e3 * t_obj / reps, 1e3 * t_dyn / reps); + + // End-to-end L-BFGS relaxation shaped like a qew_static run + // (preconditioned, flat start). Tracks the minimiser-loop overheads + // (line-search objective evals, two-loop recursion, history updates) + // that the per-kernel timings above cannot see. + const int snx = 4096, sny = 512; + const real_t sLx = 32.0, sLy = 1.0, sxi = 0.1; + FilteredNoise snoise(snx, sny, sLx, sLy, 1.0, sxi, sxi, /*seed=*/1); + Params sp; + sp.physical_size = sLx; + sp.line_tension = 1.0; // lambda_p/xi = 1 + sp.driving_force = 0.0; + sp.model = Model::Arclength; + View1D hs("hs", snx); + Kokkos::deep_copy(hs, static_cast(0.5) * sLy); + LbfgsParams sopt; + sopt.ftol = 1e-6; + sopt.max_iter = 5000; + sopt.precondition = true; + sopt.precond_shift = + static_cast(0.5) * (sLx / snx) / (sxi * sxi); + timer.reset(); + const LbfgsResult sres = + lbfgs_minimize(hs, sp, snoise.device_noise(), sopt); + const double t_lbfgs = timer.seconds(); + std::printf( + "lbfgs(n=%d, precond): %6.3fs iters=%d feval=%ld %s\n", snx, + t_lbfgs, sres.iterations, sres.n_func_evals, + sres.converged ? "converged" : "NOT converged"); } Kokkos::finalize(); return 0; diff --git a/executables/qew_depinning/main.cpp b/executables/qew_depinning/main.cpp index b171b5d..2242056 100644 --- a/executables/qew_depinning/main.cpp +++ b/executables/qew_depinning/main.cpp @@ -31,8 +31,8 @@ int main(int argc, char *argv[]) { cfg.model = Model::Arclength; const char *model_name = "arclength"; cfg.amplitude = 1.0; - cfg.xi = 0.1; // correlation length - cfg.pinning_length = 5.0; // lambda_p/xi (stiff) + cfg.xi = 0.1; // correlation length + cfg.pinning_length_over_xi = 5.0; // lambda_p/xi (stiff) cfg.Lx = 8.0; cfg.Ly = 1.0; cfg.nx = 1024; @@ -52,7 +52,7 @@ int main(int argc, char *argv[]) { const DepinningResult res = depinning_ramp(cfg, noise.device_noise()); std::printf("model=%s lambda_p/xi=%.3g line_tension=%.4g\n", model_name, - static_cast(cfg.pinning_length), + static_cast(cfg.pinning_length_over_xi), static_cast(res.line_tension)); for (const DepinningStep &st : res.steps) { std::printf("f=%7.3f sweeps=%-7d mean_h=%8.4f %s\n", diff --git a/executables/qew_dynamic/main.cpp b/executables/qew_dynamic/main.cpp index 10311e2..4350c6d 100644 --- a/executables/qew_dynamic/main.cpp +++ b/executables/qew_dynamic/main.cpp @@ -30,7 +30,9 @@ int main(int argc, char *argv[]) { const char *model_name = "arclength"; cfg.amplitude = 1.0; cfg.xi = 0.1; // correlation length - cfg.pinning_length = 0.01; + // lambda_p/xi = 0.1 (floppy); previously stored as the physical length + // 0.01 and divided by xi -- the line tension is unchanged. + cfg.pinning_length_over_xi = 0.1; cfg.driving_force = 0.0; cfg.Lx = 2.0; cfg.Ly = 1.0; diff --git a/executables/qew_static/main.cpp b/executables/qew_static/main.cpp index 963b706..159a6de 100644 --- a/executables/qew_static/main.cpp +++ b/executables/qew_static/main.cpp @@ -36,7 +36,8 @@ int main(int argc, char *argv[]) { cfg.nx = 4096; cfg.ny = 512; cfg.driving_force = 0.0; - cfg.pinning_lengths = {0.01, 0.1, 1.0, 10.0, 100.0}; + // lambda_p/xi: floppy (0.01) through stiff (100), cf. the paper's sweep. + cfg.pinning_lengths_over_xi = {0.01, 0.1, 1.0, 10.0, 100.0}; cfg.ftol = 1e-6; // ample for roughness analysis cfg.max_iter = 5000; const std::uint64_t seed = 1; @@ -47,8 +48,9 @@ int main(int argc, char *argv[]) { static_sweep(cfg, noise.device_noise()); for (const StaticProfile &prof : profiles) { - std::printf("lp=%-6.3g lt=%-10.4g iters=%-7d max|F|=%.3e %s\n", - prof.pinning_length, static_cast(prof.line_tension), + std::printf("lp/xi=%-6.3g lt=%-10.4g iters=%-7d max|F|=%.3e %s\n", + prof.pinning_length_over_xi, + static_cast(prof.line_tension), prof.result.iterations, static_cast(prof.result.max_force), prof.result.converged ? "converged" : "NOT converged"); diff --git a/src/dynamics.h b/src/dynamics.h index a4030f9..fde88b9 100644 --- a/src/dynamics.h +++ b/src/dynamics.h @@ -34,27 +34,38 @@ class LangevinDynamics { template void step(const View1D &h, const Params &p, const Noise &noise, real_t drift_coeff, real_t diff_coeff) { - gradient(h, p, noise, grad_); + gradient(h, p, noise, grad_, ws_); const View1D grad = grad_; Pool pool = pool_; const real_t a = drift_coeff; const real_t b = diff_coeff; - const bool thermal = (b > 0); - Kokkos::parallel_for( - "langevin_step", static_cast(h.extent(0)), - KOKKOS_LAMBDA(const int i) { - real_t xi = 0; - if (thermal) { + const int n = static_cast(h.extent(0)); + if (b > 0) { + // One RNG state per lane, each covering a contiguous block: a + // pool acquire per chunk instead of per site (the per-site + // acquire is lock traffic that serialises on wide GPUs). + const int n_lanes = (n < 8192) ? n : 8192; + const int chunk = (n + n_lanes - 1) / n_lanes; + Kokkos::parallel_for( + "langevin_step", n_lanes, KOKKOS_LAMBDA(const int t) { auto gen = pool.get_state(); - xi = static_cast(gen.normal()); + const int lo = t * chunk; + const int hi = (lo + chunk < n) ? lo + chunk : n; + for (int i = lo; i < hi; ++i) + h(i) += -a * grad(i) + + b * static_cast(gen.normal()); pool.free_state(gen); - } - h(i) += -a * grad(i) + b * xi; - }); + }); + } else { // zero temperature: pure gradient descent, no RNG at all + Kokkos::parallel_for( + "langevin_step", n, + KOKKOS_LAMBDA(const int i) { h(i) += -a * grad(i); }); + } } private: View1D grad_; + GradientWorkspace ws_; Pool pool_; }; diff --git a/src/filtered_noise.h b/src/filtered_noise.h index 5106b1c..731cb02 100644 --- a/src/filtered_noise.h +++ b/src/filtered_noise.h @@ -29,10 +29,9 @@ namespace qew { // Cubic B-spline interpolation weights w[0..3] for the four neighbours -// {i-1, i, i+1, i+2} at fractional position t in [0,1), and their derivatives -// dw[0..3] with respect to t (i.e. d/d(grid coordinate)). +// {i-1, i, i+1, i+2} at fractional position t in [0,1). KOKKOS_INLINE_FUNCTION -void bspline_weights(real_t t, real_t w[4], real_t dw[4]) { +void bspline_weights_value(real_t t, real_t w[4]) { const real_t t2 = t * t; const real_t t3 = t2 * t; const real_t omt = static_cast(1) - t; @@ -44,7 +43,16 @@ void bspline_weights(real_t t, real_t w[4], real_t dw[4]) { static_cast(3) * t2 - static_cast(3) * t3) / static_cast(6); w[3] = t3 / static_cast(6); +} + +// As above, plus the weight derivatives dw[0..3] with respect to t (i.e. +// d/d(grid coordinate)). +KOKKOS_INLINE_FUNCTION +void bspline_weights(real_t t, real_t w[4], real_t dw[4]) { + bspline_weights_value(t, w); + const real_t t2 = t * t; + const real_t omt = static_cast(1) - t; dw[0] = -omt * omt / static_cast(2); dw[1] = (-static_cast(12) * t + static_cast(9) * t2) / static_cast(6); @@ -63,21 +71,40 @@ struct DeviceNoise { real_t dx = 1; // physical grid spacing in x real_t dy = 1; // physical grid spacing in y + // Physical -> grid coordinates: wrap into [0, n), split into the base cell + // index and fraction, and build the four periodic neighbour indices with + // conditional wraps (integer % is expensive in these hot device loops). KOKKOS_INLINE_FUNCTION - NoiseSample sample(real_t x, real_t y) const { - // Physical -> grid coordinates, wrapped into [0, n) (positive modulo). + void locate(real_t x, real_t y, int ixn[4], int iyn[4], real_t &fx, + real_t &fy) const { real_t gx = x / dx; real_t gy = y / dy; gx -= Kokkos::floor(gx / static_cast(nx)) * static_cast(nx); gy -= Kokkos::floor(gy / static_cast(ny)) * static_cast(ny); - const int ix = static_cast(Kokkos::floor(gx)); - const int iy = static_cast(Kokkos::floor(gy)); - const real_t fx = gx - static_cast(ix); - const real_t fy = gy - static_cast(iy); + int ix = static_cast(Kokkos::floor(gx)); + int iy = static_cast(Kokkos::floor(gy)); + // The wrap can round to exactly n; fold it back before the fraction. + if (ix >= nx) { ix -= nx; gx -= static_cast(nx); } + if (iy >= ny) { iy -= ny; gy -= static_cast(ny); } + fx = gx - static_cast(ix); + fy = gy - static_cast(iy); + + ixn[0] = (ix == 0) ? nx - 1 : ix - 1; + ixn[1] = ix; + ixn[2] = (ix + 1 == nx) ? 0 : ix + 1; + ixn[3] = (ixn[2] + 1 == nx) ? 0 : ixn[2] + 1; + iyn[0] = (iy == 0) ? ny - 1 : iy - 1; + iyn[1] = iy; + iyn[2] = (iy + 1 == ny) ? 0 : iy + 1; + iyn[3] = (iyn[2] + 1 == ny) ? 0 : iyn[2] + 1; + } - const int ixn[4] = {(ix - 1 + nx) % nx, ix % nx, (ix + 1) % nx, (ix + 2) % nx}; - const int iyn[4] = {(iy - 1 + ny) % ny, iy % ny, (iy + 1) % ny, (iy + 2) % ny}; + KOKKOS_INLINE_FUNCTION + NoiseSample sample(real_t x, real_t y) const { + int ixn[4], iyn[4]; + real_t fx, fy; + locate(x, y, ixn, iyn, fx, fy); real_t wx[4], dwx[4], wy[4], dwy[4]; bspline_weights(fx, wx, dwx); @@ -105,6 +132,48 @@ struct DeviceNoise { s.dv_dy = dv_dgy / dy; return s; } + + // Value only -- the objective() hot path (line-search trials) does not + // need the derivatives, so skip the dw weights and both accumulators. + KOKKOS_INLINE_FUNCTION + real_t sample_value(real_t x, real_t y) const { + int ixn[4], iyn[4]; + real_t fx, fy; + locate(x, y, ixn, iyn, fx, fy); + + real_t wx[4], wy[4]; + bspline_weights_value(fx, wx); + bspline_weights_value(fy, wy); + + real_t v = 0; + for (int a = 0; a < 4; ++a) { + real_t row = 0; + for (int b = 0; b < 4; ++b) row += coef(ixn[a], iyn[b]) * wy[b]; + v += wx[a] * row; + } + return v; + } + + // dV/dy only -- gradient() and the Rosso-Krauth site velocity use just the + // y-derivative at edge midpoints; dV/dx is never consumed anywhere. + KOKKOS_INLINE_FUNCTION + real_t sample_dy(real_t x, real_t y) const { + int ixn[4], iyn[4]; + real_t fx, fy; + locate(x, y, ixn, iyn, fx, fy); + + real_t wx[4], wy[4], dwy[4]; + bspline_weights_value(fx, wx); + bspline_weights(fy, wy, dwy); + + real_t dv_dgy = 0; + for (int a = 0; a < 4; ++a) { + real_t row_dy = 0; + for (int b = 0; b < 4; ++b) row_dy += coef(ixn[a], iyn[b]) * dwy[b]; + dv_dgy += wx[a] * row_dy; + } + return dv_dgy / dy; + } }; class FilteredNoise { diff --git a/src/fire.h b/src/fire.h index 9c23063..e5ecb0b 100644 --- a/src/fire.h +++ b/src/fire.h @@ -44,6 +44,7 @@ FireResult fire_minimize(const View1D &h, const Params &p, const Noise &noise, const int n = static_cast(h.extent(0)); View1D grad("fire_grad", n); View1D v("fire_v", n); + GradientWorkspace gws; Kokkos::deep_copy(v, static_cast(0)); real_t dt = (fp.dt_init > 0) ? fp.dt_init : static_cast(0.1) * fp.dt_max; @@ -62,7 +63,7 @@ FireResult fire_minimize(const View1D &h, const Params &p, const Noise &noise, return m; }; - gradient(h, p, noise, grad); + gradient(h, p, noise, grad, gws); FireResult res; for (int it = 0; it < fp.max_iter; ++it) { @@ -120,7 +121,7 @@ FireResult fire_minimize(const View1D &h, const Params &p, const Noise &noise, Kokkos::parallel_for( "fire_xstep", n, KOKKOS_LAMBDA(const int i) { h(i) += dt * v(i); }); - gradient(h, p, noise, grad); + gradient(h, p, noise, grad, gws); } res.iterations = fp.max_iter; diff --git a/src/lbfgs.h b/src/lbfgs.h index f5c5979..a3b2333 100644 --- a/src/lbfgs.h +++ b/src/lbfgs.h @@ -3,6 +3,7 @@ #include #include +#include #include #include @@ -101,7 +102,8 @@ LbfgsResult lbfgs_minimize(const View1D &h, const Params &p, const Noise &noise, const int m = lp.m; View1D g("lbfgs_g", n), g_new("lbfgs_gnew", n), d("lbfgs_d", n); - View1D q("lbfgs_q", n), x0("lbfgs_x0", n), s_tmp("lbfgs_s", n), y_tmp("lbfgs_y", n); + View1D q("lbfgs_q", n), s_tmp("lbfgs_s", n), y_tmp("lbfgs_y", n); + GradientWorkspace gws; std::vector S(m), Y(m); for (int k = 0; k < m; ++k) { S[k] = View1D("lbfgs_S", n); @@ -143,9 +145,11 @@ LbfgsResult lbfgs_minimize(const View1D &h, const Params &p, const Noise &noise, irfft_plan.emplace(exec, qhat, d, KokkosFFT::Direction::backward, -1); } - gradient(h, p, noise, g); + gradient(h, p, noise, g, gws); LbfgsResult res; int ls_fail_streak = 0; + real_t phi = 0; // objective at the current iterate, carried across + // iterations (each accepted trial evaluates it anyway) for (int it = 0; it < lp.max_iter; ++it) { const real_t gmax = lbfgs_max_abs(n, g); @@ -196,9 +200,13 @@ LbfgsResult lbfgs_minimize(const View1D &h, const Params &p, const Noise &noise, } // ---- Armijo backtracking line search ---- - Kokkos::deep_copy(x0, h); - const real_t phi0 = objective(h, p, noise); // h == x0 here - ++res.n_func_evals; + // Trials are evaluated at h + alpha*d without writing h; the step is + // committed only on acceptance, so failure needs no saved copy/undo. + if (it == 0) { // later iterations reuse the accepted trial value + phi = objective(h, p, noise); + ++res.n_func_evals; + } + const real_t phi0 = phi; // objective at the current iterate // alpha0 = 1 for L-BFGS / preconditioned steps (well-scaled direction); // for an unpreconditioned first steepest-descent step, scale by 1/|g|. real_t alpha = static_cast(1); @@ -208,20 +216,16 @@ LbfgsResult lbfgs_minimize(const View1D &h, const Params &p, const Noise &noise, } bool ok = false; for (int ls = 0; ls < lp.max_ls; ++ls) { - const real_t al = alpha; - Kokkos::parallel_for( - "lbfgs_trial", n, - KOKKOS_LAMBDA(const int i) { h(i) = x0(i) + al * d(i); }); - const real_t phi = objective(h, p, noise); + const real_t phi_trial = objective(h, d, alpha, p, noise); ++res.n_func_evals; - if (phi <= phi0 + lp.c1 * alpha * dphi0) { + if (phi_trial <= phi0 + lp.c1 * alpha * dphi0) { + phi = phi_trial; ok = true; break; } alpha *= lp.ls_shrink; } if (!ok) { - Kokkos::deep_copy(h, x0); // undo the failed trial if (slots.empty() || ++ls_fail_streak > 2) { res.iterations = it; res.max_force = gmax; @@ -234,17 +238,27 @@ LbfgsResult lbfgs_minimize(const View1D &h, const Params &p, const Noise &noise, } ls_fail_streak = 0; - // ---- curvature pair s = alpha d, y = g_new - g ---- - gradient(h, p, noise, g_new); + // Commit the accepted step. const real_t al = alpha; Kokkos::parallel_for( - "lbfgs_sy", n, KOKKOS_LAMBDA(const int i) { - s_tmp(i) = al * d(i); - y_tmp(i) = g_new(i) - g(i); - }); - const real_t sy = lbfgs_dot(n, s_tmp, y_tmp); - const real_t ss = lbfgs_dot(n, s_tmp, s_tmp); - const real_t yy = lbfgs_dot(n, y_tmp, y_tmp); + "lbfgs_step", n, KOKKOS_LAMBDA(const int i) { h(i) += al * d(i); }); + + // ---- curvature pair s = alpha d, y = g_new - g ---- + // One fused kernel writes the pair and reduces all three scalars. + gradient(h, p, noise, g_new, gws); + real_t sy = 0, ss = 0, yy = 0; + Kokkos::parallel_reduce( + "lbfgs_sy", n, + KOKKOS_LAMBDA(const int i, real_t &asy, real_t &ass, real_t &ayy) { + const real_t s = al * d(i); + const real_t y = g_new(i) - g(i); + s_tmp(i) = s; + y_tmp(i) = y; + asy += s * y; + ass += s * s; + ayy += y * y; + }, + sy, ss, yy); if (sy > lp.curvature_eps * std::sqrt(ss * yy)) { // skip safeguard int row; if (static_cast(slots.size()) < m) { @@ -253,15 +267,17 @@ LbfgsResult lbfgs_minimize(const View1D &h, const Params &p, const Noise &noise, row = slots.front(); slots.erase(slots.begin()); } - Kokkos::deep_copy(S[row], s_tmp); - Kokkos::deep_copy(Y[row], y_tmp); + // Swap the pair into the history slot (s_tmp/y_tmp become the + // slot's old buffers, reused as scratch next iteration). + std::swap(S[row], s_tmp); + std::swap(Y[row], y_tmp); rho[row] = static_cast(1) / sy; gamma = sy / yy; slots.push_back(row); } else { ++res.n_curvature_skips; // non-positive curvature: keep H0 pos. def. } - Kokkos::deep_copy(g, g_new); + std::swap(g, g_new); // g <- g_new; old g becomes scratch } res.iterations = lp.max_iter; diff --git a/src/qew_io.h b/src/qew_io.h index 4da78d7..bdb7f7b 100644 --- a/src/qew_io.h +++ b/src/qew_io.h @@ -19,13 +19,16 @@ namespace qew { +// NB: the "pinning_length" keys/values carry the dimensionless ratio +// lambda_p/xi, the paper's control parameter (see qew_runs.h). inline void write_static(HighFive::File &file, const StaticConfig &c, const std::vector &profiles) { file.createDataSet("physical_size", static_cast(c.Lx)); for (const StaticProfile &prof : profiles) { - const std::string prefix = "driving_force=" + fmt("%.1g", c.driving_force) + - "/pinning_length=" + fmt("%.3g", prof.pinning_length); - file.createDataSet(prefix + "/pinning_length", prof.pinning_length); + const std::string prefix = + "driving_force=" + fmt("%.1g", c.driving_force) + + "/pinning_length=" + fmt("%.3g", prof.pinning_length_over_xi); + file.createDataSet(prefix + "/pinning_length", prof.pinning_length_over_xi); file.createDataSet(prefix + "/h", prof.h); } } @@ -39,7 +42,7 @@ inline void write_dynamic(HighFive::File &file, const DynamicConfig &c, "/driving_force=" + fmt("%.1g", c.driving_force) + "/line_tension=" + fmt("%.3g", res.line_tension); file.createDataSet(prefix + "/pinning_length", - static_cast(c.pinning_length)); + static_cast(c.pinning_length_over_xi)); file.createDataSet(prefix + "/h", snap.h); } } @@ -52,7 +55,7 @@ inline void write_depinning(HighFive::File &file, const DepinningConfig &c, const std::string prefix = "driving_force=" + fmt("%.4g", st.f); file.createDataSet(prefix + "/h", st.h); file.createDataSet(prefix + "/pinning_length", - static_cast(c.pinning_length)); + static_cast(c.pinning_length_over_xi)); } if (res.f_c_upper_bound >= 0) file.createDataSet("f_c_upper_bound", diff --git a/src/qew_model.h b/src/qew_model.h index c4b0eb1..e9e7257 100644 --- a/src/qew_model.h +++ b/src/qew_model.h @@ -26,12 +26,32 @@ struct Params { Model model = Model::Linear; }; -// Total energy E[h] = sum_i dx * (line_energy_i + V_i - f h_i), where the -// elastic and noise terms are evaluated on the forward edge (i, i+1) exactly as -// in qew_objective. Runs as a parallel_reduce on the default execution space. +namespace detail { + +// Noise access points: prefer the cheap single-quantity samplers when the +// functor provides them (DeviceNoise does); analytic test noises only supply +// the full sample() and fall back to it. Keeps one call site per kernel. template -real_t objective(const View1D &h, const Params &p, const Noise &noise) { - const int n = static_cast(h.extent(0)); +KOKKOS_INLINE_FUNCTION real_t noise_value(const Noise &noise, real_t x, real_t y) { + if constexpr (requires { noise.sample_value(x, y); }) + return noise.sample_value(x, y); + else + return noise.sample(x, y).v; +} + +template +KOKKOS_INLINE_FUNCTION real_t noise_dvdy(const Noise &noise, real_t x, real_t y) { + if constexpr (requires { noise.sample_dy(x, y); }) + return noise.sample_dy(x, y); + else + return noise.sample(x, y).dv_dy; +} + +// Shared objective body: the line position enters only through the +// device-callable functor `at(i)`, so the plain evaluation and the +// line-search trial evaluation share one analytic definition. +template +real_t objective_eval(int n, const Params &p, const Noise &noise, const At &at) { const real_t dx = p.physical_size / static_cast(n); const real_t lt = p.line_tension; const real_t f = p.driving_force; @@ -41,8 +61,10 @@ real_t objective(const View1D &h, const Params &p, const Noise &noise) { Kokkos::parallel_reduce( "qew::objective", n, KOKKOS_LAMBDA(const int i, real_t &acc) { - const int ip = (i + 1) % n; // periodic forward neighbour - const real_t dh_dx = (h(ip) - h(i)) / dx; + const int ip = (i + 1 == n) ? 0 : i + 1; // periodic forward neighbour + const real_t hi = at(i); + const real_t hip = at(ip); + const real_t dh_dx = (hip - hi) / dx; real_t line_energy; if (model == Model::Linear) { @@ -55,61 +77,129 @@ real_t objective(const View1D &h, const Params &p, const Noise &noise) { // Edge midpoint, matching python xcenter/hcenter (periodic wrap). const real_t xc = static_cast(i + ip) * static_cast(0.5) * dx; - const real_t hc = (h(ip) + h(i)) * static_cast(0.5); - const NoiseSample s = noise.sample(xc, hc); - - acc += dx * (line_energy + s.v - f * h(i)); + const real_t hc = (hip + hi) * static_cast(0.5); + acc += dx * (line_energy + noise_value(noise, xc, hc) - f * hi); }, total); return total; } +struct LineAt { + View1D h; + KOKKOS_INLINE_FUNCTION real_t operator()(int i) const { return h(i); } +}; + +struct TrialAt { + View1D x0, d; + real_t alpha; + KOKKOS_INLINE_FUNCTION real_t operator()(int i) const { + return x0(i) + alpha * d(i); + } +}; + +} // namespace detail + +// Total energy E[h] = sum_i dx * (line_energy_i + V_i - f h_i), where the +// elastic and noise terms are evaluated on the forward edge (i, i+1) exactly as +// in qew_objective. Runs as a parallel_reduce on the default execution space. +template +real_t objective(const View1D &h, const Params &p, const Noise &noise) { + return detail::objective_eval(static_cast(h.extent(0)), p, noise, + detail::LineAt{h}); +} + +// E at the line-search trial point x0 + alpha*d, evaluated without +// materialising it -- saves a kernel launch and n stores per trial, and the +// caller commits (or discards) the step afterwards. +template +real_t objective(const View1D &x0, const View1D &d, real_t alpha, + const Params &p, const Noise &noise) { + return detail::objective_eval(static_cast(x0.extent(0)), p, noise, + detail::TrialAt{x0, d, alpha}); +} + +// Scratch for the two-pass gradient: per-edge quantities computed once and +// shared by the two adjacent sites (the dominant bicubic gather would +// otherwise run twice per edge). Own one of these in loops that call +// gradient() repeatedly so the allocations are reused; the workspace-free +// overload below allocates on the fly. +struct GradientWorkspace { + View1D edge_dvdy; // dV/dy at the midpoint of edge (i, i+1) + View1D edge_slope; // arclength only: dh/sqrt(1+dh^2) on edge (i, i+1) + + void ensure(int n, Model model) { + if (static_cast(edge_dvdy.extent(0)) != n) + edge_dvdy = View1D("qew_edge_dvdy", n); + if (model == Model::Arclength && + static_cast(edge_slope.extent(0)) != n) + edge_slope = View1D("qew_edge_slope", n); + } +}; + // Force grad_i = dE/dh_i. Elastic term: harmonic Laplacian (linear) or // slope-normalized difference (arclength); noise term: half the sum of dV/dy at -// the left edge (i-1, i) and right edge (i, i+1). Exact port of qew_gradient, -// run as a parallel_for. `grad` must be pre-allocated with the same extent as h. +// the left edge (i-1, i) and right edge (i, i+1). Exact port of qew_gradient. +// `grad` must be pre-allocated with the same extent as h. +// +// Two passes: edge kernel (one noise gather + one slope per edge), then a +// cheap per-site combine. Site i's left edge is edge i-1, its right edge is +// edge i; the edge midpoints match python xcenter/hcenter (periodic wrap), so +// the result is identical to the former one-pass version that sampled every +// edge twice. template void gradient(const View1D &h, const Params &p, const Noise &noise, - const View1D &grad) { + const View1D &grad, GradientWorkspace &ws) { const int n = static_cast(h.extent(0)); const real_t dx = p.physical_size / static_cast(n); const real_t lt = p.line_tension; const real_t f = p.driving_force; const Model model = p.model; + ws.ensure(n, model); + const View1D edvdy = ws.edge_dvdy; + const View1D eslope = ws.edge_slope; + const bool arc = (model == Model::Arclength); + + Kokkos::parallel_for( + "qew::gradient_edges", n, KOKKOS_LAMBDA(const int i) { + const int ip = (i + 1 == n) ? 0 : i + 1; + const real_t xc = + static_cast(i + ip) * static_cast(0.5) * dx; + const real_t hc = (h(ip) + h(i)) * static_cast(0.5); + edvdy(i) = detail::noise_dvdy(noise, xc, hc); + if (arc) { + const real_t dh = (h(ip) - h(i)) / dx; + eslope(i) = + dh / Kokkos::sqrt(static_cast(1) + dh * dh); + } + }); + Kokkos::parallel_for( "qew::gradient", n, KOKKOS_LAMBDA(const int i) { - const int ip = (i + 1) % n; // right neighbour - const int im = (i - 1 + n) % n; // left neighbour + const int ip = (i + 1 == n) ? 0 : i + 1; // right neighbour + const int im = (i == 0) ? n - 1 : i - 1; // left neighbour real_t line_force; if (model == Model::Linear) { const real_t d2h = h(im) + h(ip) - static_cast(2) * h(i); line_force = -lt * d2h / (dx * dx); } else { - const real_t dhr = (h(ip) - h(i)) / dx; // right slope - const real_t dhl = (h(i) - h(im)) / dx; // left slope - line_force = - lt * - (dhl / Kokkos::sqrt(static_cast(1) + dhl * dhl) - - dhr / Kokkos::sqrt(static_cast(1) + dhr * dhr)) / - dx; + line_force = lt * (eslope(im) - eslope(i)) / dx; } - // dV/dy at the left and right edge midpoints (periodic wrap matches - // python xcenter_left/right and hcenter_left/right). - const real_t xcl = static_cast(im + i) * static_cast(0.5) * dx; - const real_t hcl = (h(im) + h(i)) * static_cast(0.5); - const real_t xcr = static_cast(i + ip) * static_cast(0.5) * dx; - const real_t hcr = (h(ip) + h(i)) * static_cast(0.5); - const NoiseSample sl = noise.sample(xcl, hcl); - const NoiseSample sr = noise.sample(xcr, hcr); - grad(i) = dx * (line_force + - static_cast(0.5) * (sl.dv_dy + sr.dv_dy) - f); + static_cast(0.5) * + (edvdy(im) + edvdy(i)) - f); }); } +template +void gradient(const View1D &h, const Params &p, const Noise &noise, + const View1D &grad) { + GradientWorkspace ws; + gradient(h, p, noise, grad, ws); +} + } // namespace qew #endif // QEW_MODEL_H diff --git a/src/qew_runs.h b/src/qew_runs.h index 486cc25..1bc5d18 100644 --- a/src/qew_runs.h +++ b/src/qew_runs.h @@ -22,7 +22,16 @@ // the kernels it drives, so tests can pass a cheap analytic potential instead of // paying for a FilteredNoise FFT setup. The loop bodies are a faithful copy of // the original mains -- behaviour (timestep heuristic, anneal schedule, force -// ramp + f_c bracketing, pinning-length -> line-tension map) is unchanged. +// ramp + f_c bracketing) is unchanged. +// +// Pinning-length convention (Le et al., arXiv:2410.21838, eq. (10)): the Larkin +// length obeys lambda_p / xi = (Gamma / U0)^{2/3}, varied by holding the +// disorder (amplitude U0, correlation length xi) fixed and scanning the line +// tension Gamma. All drivers therefore take the DIMENSIONLESS ratio +// lambda_p/xi and invert it as +// Gamma = amplitude * (lambda_p/xi)^{3/2}. +// The ratio is also what the HDF5 "pinning_length" keys/values carry, matching +// the paper's labelling (lambda_p/xi = 0.1 floppy ... 100 stiff). namespace qew { inline std::vector to_host_vec(const View1D &h) { @@ -39,19 +48,20 @@ inline std::vector to_host_vec(const View1D &h) { struct StaticConfig { Model model = Model::Arclength; real_t amplitude = 1.0; - real_t xi = 0.1; // correlation length + real_t xi = 0.1; // correlation length; MUST match the noise's xi real_t Lx = 32.0; real_t Ly = 1.0; int nx = 4096; int ny = 512; // noise y-resolution (used only at noise setup) real_t driving_force = 0.0; - std::vector pinning_lengths = {0.01, 0.1, 1.0, 10.0, 100.0}; + // lambda_p/xi ratios to sweep, floppy (<1) to stiff (>1). + std::vector pinning_lengths_over_xi = {0.01, 0.1, 1.0, 10.0, 100.0}; real_t ftol = 1e-6; int max_iter = 5000; }; struct StaticProfile { - double pinning_length = 0; + double pinning_length_over_xi = 0; real_t line_tension = 0; std::vector h; LbfgsResult result; @@ -60,8 +70,8 @@ struct StaticProfile { template std::vector static_sweep(const StaticConfig &c, const Noise &noise) { std::vector out; - out.reserve(c.pinning_lengths.size()); - for (double lp : c.pinning_lengths) { + out.reserve(c.pinning_lengths_over_xi.size()); + for (double lp : c.pinning_lengths_over_xi) { Params p; p.physical_size = c.Lx; p.line_tension = std::pow(lp, static_cast(1.5)) * c.amplitude; @@ -82,7 +92,7 @@ std::vector static_sweep(const StaticConfig &c, const Noise &nois const LbfgsResult r = lbfgs_minimize(h, p, noise, opt); StaticProfile prof; - prof.pinning_length = lp; + prof.pinning_length_over_xi = lp; prof.line_tension = p.line_tension; prof.h = to_host_vec(h); prof.result = r; @@ -97,8 +107,8 @@ std::vector static_sweep(const StaticConfig &c, const Noise &nois struct DynamicConfig { Model model = Model::Arclength; real_t amplitude = 1.0; - real_t xi = 0.1; - real_t pinning_length = 0.01; + real_t xi = 0.1; // correlation length; MUST match the noise's xi + real_t pinning_length_over_xi = 0.1; // lambda_p/xi (floppy) real_t driving_force = 0.0; real_t Lx = 2.0; real_t Ly = 1.0; @@ -131,7 +141,7 @@ DynamicResult dynamic_anneal(const DynamicConfig &c, const Noise &noise) { const real_t dt = std::min(static_cast(0.1), tau / 10 * std::min(fac, fac * fac)); const real_t line_tension = - std::pow(c.pinning_length / c.xi, static_cast(1.5)) * c.amplitude; + std::pow(c.pinning_length_over_xi, static_cast(1.5)) * c.amplitude; Params p; p.physical_size = c.Lx; @@ -172,8 +182,8 @@ DynamicResult dynamic_anneal(const DynamicConfig &c, const Noise &noise) { struct DepinningConfig { Model model = Model::Arclength; real_t amplitude = 1.0; - real_t xi = 0.1; - real_t pinning_length = 5.0; + real_t xi = 0.1; // correlation length; MUST match the noise's xi + real_t pinning_length_over_xi = 5.0; // lambda_p/xi (stiff) real_t Lx = 8.0; real_t Ly = 1.0; int nx = 1024; @@ -198,7 +208,7 @@ struct DepinningResult { template DepinningResult depinning_ramp(const DepinningConfig &c, const Noise &noise) { const real_t line_tension = - std::pow(c.pinning_length, static_cast(1.5)) * c.amplitude; + std::pow(c.pinning_length_over_xi, static_cast(1.5)) * c.amplitude; Params p; p.physical_size = c.Lx; @@ -211,7 +221,10 @@ DepinningResult depinning_ramp(const DepinningConfig &c, const Noise &noise) { DepinningResult res; res.line_tension = line_tension; - for (real_t f = 0; f <= c.f_max; f += c.f_step) { + // Integer counter: accumulating f += f_step drifts for non-representable + // steps and can drop the final force (and perturb the HDF5 key strings). + for (int k = 0; static_cast(k) * c.f_step <= c.f_max; ++k) { + const real_t f = static_cast(k) * c.f_step; p.driving_force = f; const RkResult r = rk_block(h, p, noise, c.rp); // restart from previous diff --git a/src/rosso_krauth.h b/src/rosso_krauth.h index 800ad24..87d88c1 100644 --- a/src/rosso_krauth.h +++ b/src/rosso_krauth.h @@ -1,6 +1,8 @@ #ifndef QEW_ROSSO_KRAUTH_H #define QEW_ROSSO_KRAUTH_H +#include + #include #include "qew_model.h" @@ -52,8 +54,8 @@ KOKKOS_INLINE_FUNCTION real_t rk_site_velocity(int i, int n, real_t dx, real_t l real_t f, Model model, real_t hi, real_t hl, real_t hr, const Noise &noise) { - const int im = (i - 1 + n) % n; - const int ip = (i + 1) % n; + const int im = (i == 0) ? n - 1 : i - 1; + const int ip = (i + 1 == n) ? 0 : i + 1; real_t line_force; if (model == Model::Linear) { @@ -71,8 +73,8 @@ KOKKOS_INLINE_FUNCTION real_t rk_site_velocity(int i, int n, real_t dx, real_t l const real_t hcl = (hl + hi) * static_cast(0.5); const real_t xcr = static_cast(i + ip) * static_cast(0.5) * dx; const real_t hcr = (hi + hr) * static_cast(0.5); - const real_t dvl = noise.sample(xcl, hcl).dv_dy; - const real_t dvr = noise.sample(xcr, hcr).dv_dy; + const real_t dvl = detail::noise_dvdy(noise, xcl, hcl); + const real_t dvr = detail::noise_dvdy(noise, xcr, hcr); const real_t grad = dx * (line_force + static_cast(0.5) * (dvl + dvr) - f); @@ -127,6 +129,11 @@ template RkResult rk_block(const View1D &h, const Params &p, const Noise &noise, const RkParams &rp) { const int n = static_cast(h.extent(0)); + // On a periodic ring with odd n, sites 0 and n-1 share parity but are + // nearest neighbours, so one colour pass would update both concurrently. + if (n % 2 != 0) + throw std::invalid_argument( + "rk_block: line length must be even for the red-black sweep"); const real_t dx = p.physical_size / static_cast(n); const real_t lt = p.line_tension; const real_t f = p.driving_force; @@ -141,13 +148,14 @@ RkResult rk_block(const View1D &h, const Params &p, const Noise &noise, mean_h0 /= static_cast(n); RkResult res; + real_t cum_adv = 0; // forward-only moves: mean_h = mean_h0 + cum_adv / n for (int sweep = 0; sweep < rp.max_sweeps; ++sweep) { for (int color = 0; color < 2; ++color) { Kokkos::parallel_for( "rk_sweep", n, KOKKOS_LAMBDA(const int i) { if ((i & 1) != color) return; - const int im = (i - 1 + n) % n; - const int ip = (i + 1) % n; + const int im = (i == 0) ? n - 1 : i - 1; + const int ip = (i + 1 == n) ? 0 : i + 1; bool found = false; const real_t hnew = rk_forward_zero(i, n, dx, lt, f, model, h(i), h(im), @@ -161,17 +169,14 @@ RkResult rk_block(const View1D &h, const Params &p, const Noise &noise, Kokkos::parallel_reduce( "rk_adv_sum", n, KOKKOS_LAMBDA(const int i, real_t &acc) { acc += adv(i); }, total_adv); - real_t mean_h = 0; - Kokkos::parallel_reduce( - "rk_mean", n, - KOKKOS_LAMBDA(const int i, real_t &acc) { acc += h(i); }, mean_h); - mean_h /= static_cast(n); + cum_adv += total_adv; + const real_t mean_adv = cum_adv / static_cast(n); res.sweeps = sweep + 1; res.total_advance = total_adv; - res.mean_h = mean_h; + res.mean_h = mean_h0 + mean_adv; - if (mean_h - mean_h0 > rp.runaway) { // depinned: advanced too far + if (mean_adv > rp.runaway) { // depinned: advanced too far res.runaway = true; res.blocked = false; return res; diff --git a/tests/test_executables.cpp b/tests/test_executables.cpp index 69af939..556b267 100644 --- a/tests/test_executables.cpp +++ b/tests/test_executables.cpp @@ -65,19 +65,21 @@ TEST(Executables, StaticSweepShapeAndMapping) { cfg.Ly = 1.0; cfg.nx = 16; cfg.ny = 8; - cfg.pinning_lengths = {0.1, 1.0, 4.0}; + cfg.pinning_lengths_over_xi = {0.1, 1.0, 4.0}; cfg.ftol = 1e-8; cfg.max_iter = 2000; const auto profiles = static_sweep(cfg, ZeroNoise{}); - ASSERT_EQ(profiles.size(), cfg.pinning_lengths.size()); + ASSERT_EQ(profiles.size(), cfg.pinning_lengths_over_xi.size()); for (std::size_t i = 0; i < profiles.size(); ++i) { const auto &prof = profiles[i]; - EXPECT_EQ(prof.pinning_length, cfg.pinning_lengths[i]); - // line_tension = pinning_length^1.5 * amplitude (the documented map). + EXPECT_EQ(prof.pinning_length_over_xi, cfg.pinning_lengths_over_xi[i]); + // Larkin inversion Gamma = amplitude * (lambda_p/xi)^{3/2} + // (arXiv:2410.21838 eq. (10)); shared by all three drivers. EXPECT_NEAR(prof.line_tension, - std::pow(cfg.pinning_lengths[i], 1.5) * cfg.amplitude, 1e-12); + std::pow(cfg.pinning_lengths_over_xi[i], 1.5) * cfg.amplitude, + 1e-12); ASSERT_EQ(prof.h.size(), static_cast(cfg.nx)); // Zero noise, zero force: the flat start is already the minimum. EXPECT_TRUE(prof.result.converged); @@ -94,7 +96,7 @@ TEST(Executables, DynamicAnnealTimestepAndSnapshots) { cfg.model = Model::Linear; cfg.amplitude = 1.0; cfg.xi = 0.1; - cfg.pinning_length = 0.01; + cfg.pinning_length_over_xi = 0.1; cfg.Lx = 1.0; cfg.Ly = 1.0; cfg.nx = 32; @@ -111,8 +113,9 @@ TEST(Executables, DynamicAnnealTimestepAndSnapshots) { const real_t dt_expected = std::min(static_cast(0.1), cfg.tau / 10 * std::min(fac, fac * fac)); EXPECT_NEAR(res.dt, dt_expected, 1e-12); + // Same Larkin map as the static driver: Gamma = A * (lambda_p/xi)^{3/2}. EXPECT_NEAR(res.line_tension, - std::pow(cfg.pinning_length / cfg.xi, 1.5) * cfg.amplitude, 1e-12); + std::pow(cfg.pinning_length_over_xi, 1.5) * cfg.amplitude, 1e-12); // Snapshots are logarithmically spaced: first near 10*dt, each ~10x the last. ASSERT_GE(res.snapshots.size(), 2u); @@ -137,7 +140,7 @@ TEST(Executables, DepinningRampBracketsCriticalForce) { DepinningConfig cfg; cfg.model = Model::Linear; cfg.amplitude = 1.0; - cfg.pinning_length = 1.0; + cfg.pinning_length_over_xi = 1.0; cfg.Lx = 8.0; cfg.Ly = 1.0; cfg.nx = 8; @@ -153,6 +156,10 @@ TEST(Executables, DepinningRampBracketsCriticalForce) { const DepinningResult res = depinning_ramp(cfg, noise); + // Same Larkin map as the other drivers: Gamma = A * (lambda_p/xi)^{3/2}. + EXPECT_NEAR(res.line_tension, + std::pow(cfg.pinning_length_over_xi, 1.5) * cfg.amplitude, 1e-12); + ASSERT_FALSE(res.steps.empty()); EXPECT_TRUE(res.steps.front().result.blocked); // f = 0 blocks EXPECT_EQ(res.steps.front().f, 0.0); @@ -195,7 +202,7 @@ TEST(Executables, WriteStaticKeysAndValues) { cfg.driving_force = 0.0; std::vector profiles(1); - profiles[0].pinning_length = 0.1; + profiles[0].pinning_length_over_xi = 0.1; profiles[0].line_tension = 0.5; profiles[0].h = {1.0, 2.0, 3.0, 4.0}; @@ -225,7 +232,7 @@ TEST(Executables, WriteDynamicKeysAndValues) { DynamicConfig cfg; cfg.Lx = 2.0; cfg.driving_force = 0.0; - cfg.pinning_length = 0.01; + cfg.pinning_length_over_xi = 0.1; DynamicResult res; res.line_tension = 2.0; @@ -249,7 +256,7 @@ TEST(Executables, WriteDynamicKeysAndValues) { TEST(Executables, WriteDepinningSkipsRunawayAndWritesFc) { DepinningConfig cfg; cfg.Lx = 8.0; - cfg.pinning_length = 5.0; + cfg.pinning_length_over_xi = 5.0; DepinningResult res; res.f_c_upper_bound = 0.8; diff --git a/tests/test_rosso_krauth.cpp b/tests/test_rosso_krauth.cpp index 8d2c828..bbf708b 100644 --- a/tests/test_rosso_krauth.cpp +++ b/tests/test_rosso_krauth.cpp @@ -3,6 +3,7 @@ #include #include #include +#include #include #include @@ -288,3 +289,25 @@ TEST(RossoKrauth, BlocksArclengthModel) { Kokkos::Max(max_v)); EXPECT_LT(max_v, 1e-6); } + +// An odd line length must be rejected up front: with periodic boundaries sites +// 0 and n-1 share red-black parity but are nearest neighbours, so one colour +// pass would update both concurrently (a data race, not just wrong physics). +TEST(RossoKrauth, RejectsOddLineLength) { + const int n = 17; + const SinusoidNoise noise; + + Params p; + p.physical_size = n; + p.line_tension = 1.0; + p.driving_force = 0.0; + p.model = Model::Linear; + + RkParams rp; + rp.dstep = 0.01; + rp.max_advance = 1.0; + rp.runaway = 3.0; + + View1D h = flat_line(n, 0.0); + EXPECT_THROW(rk_block(h, p, noise, rp), std::invalid_argument); +}