diff --git a/executables/qew_depinning/main.cpp b/executables/qew_depinning/main.cpp index c28d654..b171b5d 100644 --- a/executables/qew_depinning/main.cpp +++ b/executables/qew_depinning/main.cpp @@ -5,108 +5,70 @@ // stops when no blocked state exists (the line runs away) -- that brackets f_c. // Profiles are written to depinning_{model}.h5; analysis stays in Python. // -// Configuration is by the constants below, matching the Python convention. +// The ramp lives in qew::depinning_ramp (src/qew_runs.h) and the HDF5 layout in +// qew::write_depinning (src/qew_io.h); both are unit-tested. This main is just +// configuration + glue, matching the Python module-constant convention. -#include #include #include #include #include #include -#include #include "filtered_noise.h" -#include "qew_model.h" +#include "qew_io.h" +#include "qew_runs.h" #include "qew_types.h" -#include "rosso_krauth.h" using namespace qew; -namespace { -std::string fmt(const char *spec, double v) { - char buf[64]; - std::snprintf(buf, sizeof(buf), spec, v); - return std::string(buf); -} -} // namespace - int main(int argc, char *argv[]) { Kokkos::initialize(argc, argv); { // --- configuration (cf. qew_depinning.py) --- - const Model model = Model::Arclength; + DepinningConfig cfg; + cfg.model = Model::Arclength; const char *model_name = "arclength"; - const real_t amplitude = 1.0; - const real_t xi = 0.1; // correlation length - const real_t pinning_length = 5.0; // lambda_p/xi (stiff) - const real_t Lx = 8.0, Ly = 1.0; - const int nx = 1024, ny = 256; + cfg.amplitude = 1.0; + cfg.xi = 0.1; // correlation length + cfg.pinning_length = 5.0; // lambda_p/xi (stiff) + cfg.Lx = 8.0; + cfg.Ly = 1.0; + cfg.nx = 1024; + cfg.ny = 256; + cfg.f_step = 1.0; // force ramp: from 0 upward until the line depins + cfg.f_max = 40.0; + cfg.rp.dstep = static_cast(0.5) * (cfg.Ly / cfg.ny); + cfg.rp.max_advance = cfg.Ly; + cfg.rp.runaway = 2 * cfg.Ly; // advanced two periods without blocking + cfg.rp.root_tol = 1e-10; + cfg.rp.sweep_tol = 1e-8; + cfg.rp.max_sweeps = 200000; const std::uint64_t seed = 1; - const real_t line_tension = - std::pow(pinning_length, static_cast(1.5)) * amplitude; - - // Force ramp: from 0 upward until the line depins. - const real_t f_step = 1.0; - const real_t f_max = 40.0; - - FilteredNoise noise(nx, ny, Lx, Ly, amplitude, xi, xi, seed); - const DeviceNoise dn = noise.device_noise(); - - Params p; - p.physical_size = Lx; - p.line_tension = line_tension; - p.model = model; - - RkParams rp; - rp.dstep = static_cast(0.5) * (Ly / ny); - rp.max_advance = Ly; - rp.runaway = 2 * Ly; // advanced two periods without blocking = depinned - rp.root_tol = 1e-10; - rp.sweep_tol = 1e-8; - rp.max_sweeps = 200000; - - const std::string fname = std::string("depinning_") + model_name + ".h5"; - HighFive::File file(fname, HighFive::File::Truncate); - file.createDataSet("physical_size", static_cast(Lx)); + FilteredNoise noise(cfg.nx, cfg.ny, cfg.Lx, cfg.Ly, cfg.amplitude, cfg.xi, + cfg.xi, seed); + 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(pinning_length), - static_cast(line_tension)); - - View1D h("h", nx); - Kokkos::deep_copy(h, static_cast(0.5) * Ly); // flat start - - real_t f_c = -1; - for (real_t f = 0; f <= f_max; f += f_step) { - p.driving_force = f; - const RkResult r = rk_block(h, p, dn, rp); - + static_cast(cfg.pinning_length), + static_cast(res.line_tension)); + for (const DepinningStep &st : res.steps) { std::printf("f=%7.3f sweeps=%-7d mean_h=%8.4f %s\n", - static_cast(f), r.sweeps, - static_cast(r.mean_h), - r.blocked ? "blocked" - : (r.runaway ? "DEPINNED" : "timeout")); - - if (r.blocked) { - auto hh = Kokkos::create_mirror_view(h); - Kokkos::deep_copy(hh, h); - std::vector hvec(nx); - for (int i = 0; i < nx; ++i) hvec[i] = hh(i); - - const std::string prefix = "driving_force=" + fmt("%.4g", f); - file.createDataSet(prefix + "/h", hvec); - file.createDataSet(prefix + "/pinning_length", - static_cast(pinning_length)); - } else { - f_c = f; // first force with no blocked state brackets f_c - std::printf("Line depinned near f=%.3f; stopping ramp.\n", - static_cast(f)); - break; - } + static_cast(st.f), st.result.sweeps, + static_cast(st.result.mean_h), + st.result.blocked + ? "blocked" + : (st.result.runaway ? "DEPINNED" : "timeout")); } - if (f_c >= 0) file.createDataSet("f_c_upper_bound", static_cast(f_c)); + if (res.f_c_upper_bound >= 0) + std::printf("Line depinned near f=%.3f; stopping ramp.\n", + static_cast(res.f_c_upper_bound)); + + const std::string fname = std::string("depinning_") + model_name + ".h5"; + HighFive::File file(fname, HighFive::File::Truncate); + write_depinning(file, cfg, res); std::printf("wrote %s (run tools/postprocess_static.py to add acf/psd)\n", fname.c_str()); } diff --git a/executables/qew_dynamic/main.cpp b/executables/qew_dynamic/main.cpp index 884cb12..10311e2 100644 --- a/executables/qew_dynamic/main.cpp +++ b/executables/qew_dynamic/main.cpp @@ -3,111 +3,62 @@ // annealed thermal length L_T and saves snapshots at logarithmically spaced // times to dynamic_{model}.h5. Analysis stays in Python (postprocess + plot.py). // -// Configuration is by the constants below, matching the Python convention. +// The anneal loop lives in qew::dynamic_anneal (src/qew_runs.h) and the HDF5 +// layout in qew::write_dynamic (src/qew_io.h); both are unit-tested. This main +// is just configuration + glue, matching the Python module-constant convention. -#include -#include #include #include #include #include #include -#include -#include "dynamics.h" #include "filtered_noise.h" -#include "qew_model.h" +#include "qew_io.h" +#include "qew_runs.h" #include "qew_types.h" using namespace qew; -namespace { -std::string fmt(const char *spec, double v) { - char buf[64]; - std::snprintf(buf, sizeof(buf), spec, v); - return std::string(buf); -} -} // namespace - int main(int argc, char *argv[]) { Kokkos::initialize(argc, argv); { // --- configuration (cf. qew_dynamic.py module constants) --- - const Model model = Model::Arclength; + DynamicConfig cfg; + cfg.model = Model::Arclength; const char *model_name = "arclength"; - const real_t amplitude = 1.0; - const real_t xi = 0.1; // correlation length - const real_t pinning_length = 0.01; - const real_t driving_force = 0.0; - const real_t Lx = 2.0, Ly = 1.0; - const int nx = 2048, ny = 1024; - const std::uint64_t seed = 1; - - const real_t tau = 1.0; - const real_t tl_start = 0.01; - const real_t tl_end = 0.0001; - // qew_dynamic.py timestep heuristic. - const real_t fac = xi / std::max(tl_start, tl_end); - const real_t dt = std::min(static_cast(0.1), - tau / 10 * std::min(fac, fac * fac)); + cfg.amplitude = 1.0; + cfg.xi = 0.1; // correlation length + cfg.pinning_length = 0.01; + cfg.driving_force = 0.0; + cfg.Lx = 2.0; + cfg.Ly = 1.0; + cfg.nx = 2048; + cfg.ny = 1024; + cfg.seed = 1; + cfg.tau = 1.0; + cfg.tl_start = 0.01; + cfg.tl_end = 0.0001; // Production qew_dynamic.py uses total_time = 1e6; kept short here so the // demo finishes quickly. Increase for a real anneal. - const real_t total_time = 1000.0; + cfg.total_time = 1000.0; - const real_t line_tension = - std::pow(pinning_length / xi, static_cast(1.5)) * amplitude; + FilteredNoise noise(cfg.nx, cfg.ny, cfg.Lx, cfg.Ly, cfg.amplitude, cfg.xi, + cfg.xi, cfg.seed); + const DynamicResult res = dynamic_anneal(cfg, noise.device_noise()); - FilteredNoise noise(nx, ny, Lx, Ly, amplitude, xi, xi, seed); - const DeviceNoise dn = noise.device_noise(); - - Params p; - p.physical_size = Lx; - p.line_tension = line_tension; - p.driving_force = driving_force; - p.model = model; - - View1D h("h", nx); - Kokkos::deep_copy(h, static_cast(0.5) * Ly); // flat start - LangevinDynamics dyn(nx, seed + 1); + std::printf("dt=%.3g total_time=%.3g steps=%ld line_tension=%.4g\n", + static_cast(res.dt), + static_cast(res.total_time), + static_cast(res.total_time / res.dt), + static_cast(res.line_tension)); + for (const DynamicSnapshot &snap : res.snapshots) + std::printf(" saved t=%.6g\n", static_cast(snap.t)); const std::string fname = std::string("dynamic_") + model_name + ".h5"; HighFive::File file(fname, HighFive::File::Truncate); - file.createDataSet("physical_size", static_cast(Lx)); - - std::printf("dt=%.3g total_time=%.3g steps=%ld line_tension=%.4g\n", - static_cast(dt), static_cast(total_time), - static_cast(total_time / dt), - static_cast(line_tension)); - - auto save = [&](real_t t) { - auto hh = Kokkos::create_mirror_view(h); - Kokkos::deep_copy(hh, h); - std::vector hvec(nx); - for (int i = 0; i < nx; ++i) hvec[i] = hh(i); - const std::string prefix = "timestep=" + fmt("%.6g", t) + - "/driving_force=" + fmt("%.1g", driving_force) + - "/line_tension=" + fmt("%.3g", line_tension); - file.createDataSet(prefix + "/pinning_length", pinning_length); - file.createDataSet(prefix + "/h", hvec); - std::printf(" saved t=%.6g\n", static_cast(t)); - }; - - real_t t = 0; - real_t next_save = 10 * dt; - real_t tl = tl_start; - while (t < total_time) { - const real_t drift = tl * dt / (tau * amplitude); - const real_t diff = tl * std::sqrt(2 * dt / tau); - dyn.step(h, p, dn, drift, diff); - - tl = tl_start + (tl_end - tl_start) * (t / total_time); // anneal - t += dt; - if (t > next_save - dt / 2) { - save(t); - next_save *= 10; - } - } + write_dynamic(file, cfg, res); std::printf( "wrote %s (run tools/postprocess_static.py to add acf/psd)\n", fname.c_str()); diff --git a/executables/qew_static/main.cpp b/executables/qew_static/main.cpp index ddb9633..963b706 100644 --- a/executables/qew_static/main.cpp +++ b/executables/qew_static/main.cpp @@ -1,94 +1,62 @@ -// Static qEW: FIRE energy minimisation sweeping the pinning length, the C++/GPU -// port of ../../Overleaf/qEW/qew_static.py. Relaxes a flat line to a local -// minimum at zero driving force for each pinning length and writes the profiles -// to static_{model}.h5. Analysis (ACF/PSD/Hurst) stays in Python: run +// Static qEW: L-BFGS energy minimisation sweeping the pinning length, the +// C++/GPU port of ../../Overleaf/qEW/qew_static.py. Relaxes a flat line to a +// local minimum at zero driving force for each pinning length and writes the +// profiles to static_{model}.h5. Analysis (ACF/PSD/Hurst) stays in Python: run // tools/postprocess_static.py on the output to add the datasets plot.py reads. // -// Configuration is by the constants below, matching the Python convention. +// The sweep itself lives in qew::static_sweep (src/qew_runs.h) and the HDF5 +// layout in qew::write_static (src/qew_io.h); both are unit-tested. This main is +// just the configuration + glue, matching the Python module-constant convention. -#include #include #include #include #include #include -#include #include "filtered_noise.h" -#include "lbfgs.h" -#include "qew_model.h" +#include "qew_io.h" +#include "qew_runs.h" #include "qew_types.h" using namespace qew; -namespace { -std::string fmt(const char *spec, double v) { - char buf[64]; - std::snprintf(buf, sizeof(buf), spec, v); - return std::string(buf); -} -} // namespace - int main(int argc, char *argv[]) { Kokkos::initialize(argc, argv); { // --- configuration (cf. qew_static.py module constants) --- - const Model model = Model::Arclength; + StaticConfig cfg; + cfg.model = Model::Arclength; const char *model_name = "arclength"; - const real_t amplitude = 1.0; - const real_t xi = 0.1; // correlation length - const real_t Lx = 32.0, Ly = 1.0; - const int nx = 4096, ny = 512; - const real_t driving_force = 0.0; - const std::vector pinning_lengths = {0.01, 0.1, 1.0, 10.0, 100.0}; + cfg.amplitude = 1.0; + cfg.xi = 0.1; // correlation length + cfg.Lx = 32.0; + cfg.Ly = 1.0; + cfg.nx = 4096; + cfg.ny = 512; + cfg.driving_force = 0.0; + cfg.pinning_lengths = {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; - FilteredNoise noise(nx, ny, Lx, Ly, amplitude, xi, xi, seed); - const DeviceNoise dn = noise.device_noise(); + FilteredNoise noise(cfg.nx, cfg.ny, cfg.Lx, cfg.Ly, cfg.amplitude, cfg.xi, + cfg.xi, seed); + const std::vector profiles = + 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), + prof.result.iterations, + static_cast(prof.result.max_force), + prof.result.converged ? "converged" : "NOT converged"); + } const std::string fname = std::string("static_") + model_name + ".h5"; HighFive::File file(fname, HighFive::File::Truncate); - file.createDataSet("physical_size", static_cast(Lx)); - - for (double lp : pinning_lengths) { - Params p; - p.physical_size = Lx; - p.line_tension = std::pow(lp, 1.5) * amplitude; - p.driving_force = driving_force; - p.model = model; - - View1D h("h", nx); - Kokkos::deep_copy(h, static_cast(0.5) * Ly); // flat start - - LbfgsParams opt; - opt.ftol = 1e-6; // ample for roughness analysis - opt.max_iter = 5000; - // Inverse-Laplacian FFT preconditioner: removes the stiff lines' - // ill-conditioning (cond ~ N^2). The shift is the pinning - // (non-elastic) curvature scale ~ dx*A/xi^2, so the preconditioner - // degrades to a well-scaled scalar step on floppy lines (where the - // elastic term is negligible) and is the true Laplacian inverse on - // stiff ones. - opt.precondition = true; - opt.precond_shift = static_cast(0.5) * (Lx / nx) * amplitude / (xi * xi); - const LbfgsResult r = lbfgs_minimize(h, p, dn, opt); - - auto hh = Kokkos::create_mirror_view(h); - Kokkos::deep_copy(hh, h); - std::vector hvec(nx); - for (int i = 0; i < nx; ++i) hvec[i] = hh(i); - - const std::string prefix = "driving_force=" + fmt("%.1g", driving_force) + - "/pinning_length=" + fmt("%.3g", lp); - file.createDataSet(prefix + "/pinning_length", lp); - file.createDataSet(prefix + "/h", hvec); - - std::printf("lp=%-6.3g lt=%-10.4g iters=%-7d max|F|=%.3e %s\n", lp, - static_cast(p.line_tension), r.iterations, - static_cast(r.max_force), - r.converged ? "converged" : "NOT converged"); - } + write_static(file, cfg, profiles); std::printf( "wrote %s (run tools/postprocess_static.py to add acf/psd for " "plot.py)\n", diff --git a/src/lbfgs.h b/src/lbfgs.h index 3fd965d..f5c5979 100644 --- a/src/lbfgs.h +++ b/src/lbfgs.h @@ -55,6 +55,8 @@ struct LbfgsResult { real_t max_force = 0; bool converged = false; long n_func_evals = 0; + long n_curvature_skips = 0; // (s,y) pairs rejected by the curvature safeguard + long n_ls_restarts = 0; // line-search failures recovered by history reset }; namespace detail { @@ -227,6 +229,7 @@ LbfgsResult lbfgs_minimize(const View1D &h, const Params &p, const Noise &noise, return res; } slots.clear(); // reset history, retry as steepest descent + ++res.n_ls_restarts; continue; } ls_fail_streak = 0; @@ -255,6 +258,8 @@ LbfgsResult lbfgs_minimize(const View1D &h, const Params &p, const Noise &noise, 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); } diff --git a/src/qew_format.h b/src/qew_format.h new file mode 100644 index 0000000..9df6ae6 --- /dev/null +++ b/src/qew_format.h @@ -0,0 +1,21 @@ +#ifndef QEW_FORMAT_H +#define QEW_FORMAT_H + +#include +#include + +// Tiny printf-into-std::string helper shared by the simulation executables to +// build the HDF5 group keys. The key strings (e.g. "pinning_length=%.3g") MUST +// match the Python scripts' formatting so the same analysis/plot code reads +// either output -- see test_executables.cpp, which pins the documented keys. +namespace qew { + +inline std::string fmt(const char *spec, double v) { + char buf[64]; + std::snprintf(buf, sizeof(buf), spec, v); + return std::string(buf); +} + +} // namespace qew + +#endif // QEW_FORMAT_H diff --git a/src/qew_io.h b/src/qew_io.h new file mode 100644 index 0000000..4da78d7 --- /dev/null +++ b/src/qew_io.h @@ -0,0 +1,64 @@ +#ifndef QEW_IO_H +#define QEW_IO_H + +// HDF5 output for the three simulations. Header-only (HighFive is header-only) +// and compiled only where HDF5 is available, so it is included by the +// executables and by the HDF5-guarded tests, but never by the core `lib` +// (which does not link HighFive). The dataset KEYS written here are a contract +// with the Python analysis/plot scripts -- test_executables.cpp reads them back +// and asserts the exact strings. + +#include +#include + +#include + +#include "qew_format.h" +#include "qew_runs.h" +#include "qew_types.h" + +namespace qew { + +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); + file.createDataSet(prefix + "/h", prof.h); + } +} + +inline void write_dynamic(HighFive::File &file, const DynamicConfig &c, + const DynamicResult &res) { + file.createDataSet("physical_size", static_cast(c.Lx)); + for (const DynamicSnapshot &snap : res.snapshots) { + const std::string prefix = + "timestep=" + fmt("%.6g", snap.t) + + "/driving_force=" + fmt("%.1g", c.driving_force) + + "/line_tension=" + fmt("%.3g", res.line_tension); + file.createDataSet(prefix + "/pinning_length", + static_cast(c.pinning_length)); + file.createDataSet(prefix + "/h", snap.h); + } +} + +inline void write_depinning(HighFive::File &file, const DepinningConfig &c, + const DepinningResult &res) { + file.createDataSet("physical_size", static_cast(c.Lx)); + for (const DepinningStep &st : res.steps) { + if (!st.result.blocked) continue; + 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)); + } + if (res.f_c_upper_bound >= 0) + file.createDataSet("f_c_upper_bound", + static_cast(res.f_c_upper_bound)); +} + +} // namespace qew + +#endif // QEW_IO_H diff --git a/src/qew_runs.h b/src/qew_runs.h new file mode 100644 index 0000000..486cc25 --- /dev/null +++ b/src/qew_runs.h @@ -0,0 +1,234 @@ +#ifndef QEW_RUNS_H +#define QEW_RUNS_H + +#include +#include +#include +#include +#include +#include + +#include + +#include "dynamics.h" +#include "lbfgs.h" +#include "qew_model.h" +#include "qew_types.h" +#include "rosso_krauth.h" + +// Orchestration of the three production simulations, factored out of the +// executables' main() so it can be unit-tested without HDF5 (the I/O is split +// into qew_io.h). Each routine is templated on the `Noise` functor, exactly like +// 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. +namespace qew { + +inline std::vector to_host_vec(const View1D &h) { + auto hh = Kokkos::create_mirror_view(h); + Kokkos::deep_copy(hh, h); + std::vector out(h.extent(0)); + for (std::size_t i = 0; i < out.size(); ++i) out[i] = hh(i); + return out; +} + +// --------------------------------------------------------------------------- +// Static: L-BFGS relaxation sweeping the pinning length (cf. qew_static.py). +// --------------------------------------------------------------------------- +struct StaticConfig { + Model model = Model::Arclength; + real_t amplitude = 1.0; + real_t xi = 0.1; // correlation length + 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}; + real_t ftol = 1e-6; + int max_iter = 5000; +}; + +struct StaticProfile { + double pinning_length = 0; + real_t line_tension = 0; + std::vector h; + LbfgsResult result; +}; + +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) { + Params p; + p.physical_size = c.Lx; + p.line_tension = std::pow(lp, static_cast(1.5)) * c.amplitude; + p.driving_force = c.driving_force; + p.model = c.model; + + View1D h("h", c.nx); + Kokkos::deep_copy(h, static_cast(0.5) * c.Ly); // flat start + + LbfgsParams opt; + opt.ftol = c.ftol; + opt.max_iter = c.max_iter; + // Inverse-Laplacian FFT preconditioner with the pinning-curvature shift + // ~ dx*A/xi^2 (degrades to a well-scaled scalar step on floppy lines). + opt.precondition = true; + opt.precond_shift = static_cast(0.5) * (c.Lx / c.nx) * c.amplitude / + (c.xi * c.xi); + const LbfgsResult r = lbfgs_minimize(h, p, noise, opt); + + StaticProfile prof; + prof.pinning_length = lp; + prof.line_tension = p.line_tension; + prof.h = to_host_vec(h); + prof.result = r; + out.push_back(std::move(prof)); + } + return out; +} + +// --------------------------------------------------------------------------- +// Dynamic: annealed overdamped Langevin (cf. qew_dynamic.py). +// --------------------------------------------------------------------------- +struct DynamicConfig { + Model model = Model::Arclength; + real_t amplitude = 1.0; + real_t xi = 0.1; + real_t pinning_length = 0.01; + real_t driving_force = 0.0; + real_t Lx = 2.0; + real_t Ly = 1.0; + int nx = 2048; + int ny = 1024; + std::uint64_t seed = 1; + real_t tau = 1.0; + real_t tl_start = 0.01; + real_t tl_end = 0.0001; + real_t total_time = 1000.0; +}; + +struct DynamicSnapshot { + real_t t = 0; + std::vector h; +}; + +struct DynamicResult { + real_t dt = 0; + real_t total_time = 0; + real_t line_tension = 0; + std::vector snapshots; +}; + +template +DynamicResult dynamic_anneal(const DynamicConfig &c, const Noise &noise) { + const real_t tau = c.tau; + // qew_dynamic.py timestep heuristic. + const real_t fac = c.xi / std::max(c.tl_start, c.tl_end); + 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; + + Params p; + p.physical_size = c.Lx; + p.line_tension = line_tension; + p.driving_force = c.driving_force; + p.model = c.model; + + View1D h("h", c.nx); + Kokkos::deep_copy(h, static_cast(0.5) * c.Ly); // flat start + LangevinDynamics dyn(c.nx, c.seed + 1); + + DynamicResult res; + res.dt = dt; + res.total_time = c.total_time; + res.line_tension = line_tension; + + real_t t = 0; + real_t next_save = 10 * dt; + real_t tl = c.tl_start; + while (t < c.total_time) { + const real_t drift = tl * dt / (tau * c.amplitude); + const real_t diff = tl * std::sqrt(2 * dt / tau); + dyn.step(h, p, noise, drift, diff); + + tl = c.tl_start + (c.tl_end - c.tl_start) * (t / c.total_time); // anneal + t += dt; + if (t > next_save - dt / 2) { + res.snapshots.push_back(DynamicSnapshot{t, to_host_vec(h)}); + next_save *= 10; + } + } + return res; +} + +// --------------------------------------------------------------------------- +// Depinning: Rosso-Krauth force ramp toward f_c (cf. qew_depinning.py). +// --------------------------------------------------------------------------- +struct DepinningConfig { + Model model = Model::Arclength; + real_t amplitude = 1.0; + real_t xi = 0.1; + real_t pinning_length = 5.0; + real_t Lx = 8.0; + real_t Ly = 1.0; + int nx = 1024; + int ny = 256; + real_t f_step = 1.0; + real_t f_max = 40.0; + RkParams rp; // caller fills dstep/max_advance/runaway/tolerances +}; + +struct DepinningStep { + real_t f = 0; + RkResult result; + std::vector h; // populated only when result.blocked +}; + +struct DepinningResult { + real_t line_tension = 0; + real_t f_c_upper_bound = -1; // first force with no blocked state, else -1 + std::vector steps; +}; + +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; + + Params p; + p.physical_size = c.Lx; + p.line_tension = line_tension; + p.model = c.model; + + View1D h("h", c.nx); + Kokkos::deep_copy(h, static_cast(0.5) * c.Ly); // flat start + + DepinningResult res; + res.line_tension = line_tension; + + for (real_t f = 0; f <= c.f_max; f += c.f_step) { + p.driving_force = f; + const RkResult r = rk_block(h, p, noise, c.rp); // restart from previous + + DepinningStep st; + st.f = f; + st.result = r; + if (r.blocked) { + st.h = to_host_vec(h); + } else { + res.f_c_upper_bound = f; // first unblocked force brackets f_c + } + res.steps.push_back(std::move(st)); + if (!r.blocked) break; + } + return res; +} + +} // namespace qew + +#endif // QEW_RUNS_H diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index adc92ef..6d927db 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -20,6 +20,7 @@ set(TEST_SOURCES test_dynamics.cpp test_lbfgs.cpp test_rosso_krauth.cpp + test_executables.cpp ) # Kokkos-aware GoogleTest entry point (initializes/finalizes Kokkos around the diff --git a/tests/test_dynamics.cpp b/tests/test_dynamics.cpp index 99f479c..597e6ec 100644 --- a/tests/test_dynamics.cpp +++ b/tests/test_dynamics.cpp @@ -97,3 +97,65 @@ TEST(Dynamics, EquipartitionMeanEnergy) { const double expected = 0.5 * n * kT; // (N/2) kT EXPECT_NEAR(mean_e, expected, 0.1 * expected) << "mean E = " << mean_e; } + +// Test 8c: the driving force enters the drift. At zero temperature a flat line +// in a harmonic trough at y0 plus a constant force f relaxes (elastic term +// vanishes for a flat line) to the shifted equilibrium where the site velocity +// -dE/dh = -(k(y-y0) - f) is zero, i.e. y = y0 + f/k. Exercises the f term in +// the integrated gradient (the existing tests all use f = 0). +TEST(Dynamics, DrivingForceShiftsEquilibrium) { + const int n = 32; + const real_t k = 4.0, y0 = 0.2, f = 0.6; + Params p; + p.physical_size = n; // dx = 1 + p.line_tension = 1.0; + p.driving_force = f; + p.model = Model::Linear; + HarmonicWell well{k, y0}; + + View1D h("h", n); + Kokkos::deep_copy(h, static_cast(0)); + LangevinDynamics dyn(n, /*seed=*/0); + + const real_t drift = 0.05; + for (int s = 0; s < 5000; ++s) + dyn.step(h, p, well, drift, /*diff_coeff=*/0.0); + + auto hh = Kokkos::create_mirror_view(h); + Kokkos::deep_copy(hh, h); + const real_t y_eq = y0 + f / k; + real_t max_dev = 0; + for (int i = 0; i < n; ++i) max_dev = std::max(max_dev, std::abs(hh(i) - y_eq)); + EXPECT_LT(max_dev, 1e-6); +} + +// Test 8d: with diff_coeff == 0 the RNG branch is skipped entirely, so the step +// is purely deterministic -- two integrators seeded DIFFERENTLY must produce +// identical trajectories. Pins that the thermal kick is gated on diff_coeff > 0 +// (a regression here would silently inject noise into the supposedly +// zero-temperature gradient-descent path the static cross-checks rely on). +TEST(Dynamics, ZeroDiffusionIgnoresSeed) { + const int n = 48; + Params p; + p.physical_size = n; + p.line_tension = 1.0; + p.driving_force = 0.3; + p.model = Model::Linear; + HarmonicWell well{4.0, 0.1}; + + View1D h1 = random_line(n, /*seed=*/5, -0.5, 0.5); + View1D h2 = random_line(n, /*seed=*/5, -0.5, 0.5); // same start + LangevinDynamics d1(n, /*seed=*/111); + LangevinDynamics d2(n, /*seed=*/999); // different RNG seed + + for (int s = 0; s < 50; ++s) { + d1.step(h1, p, well, 0.05, /*diff_coeff=*/0.0); + d2.step(h2, p, well, 0.05, /*diff_coeff=*/0.0); + } + + auto a = Kokkos::create_mirror_view(h1); + auto b = Kokkos::create_mirror_view(h2); + Kokkos::deep_copy(a, h1); + Kokkos::deep_copy(b, h2); + for (int i = 0; i < n; ++i) EXPECT_EQ(a(i), b(i)); +} diff --git a/tests/test_executables.cpp b/tests/test_executables.cpp new file mode 100644 index 0000000..69af939 --- /dev/null +++ b/tests/test_executables.cpp @@ -0,0 +1,283 @@ +#include "qew_runs.h" + +#include +#include +#include +#include + +#include +#include + +#include "qew_format.h" +#include "qew_types.h" + +#ifdef QEW_WITH_HDF5 +#include + +#include "qew_io.h" +#endif + +// Coverage for the simulation executables, whose orchestration was factored into +// qew::static_sweep / dynamic_anneal / depinning_ramp (src/qew_runs.h) and whose +// HDF5 layout lives in qew::write_* (src/qew_io.h). The mains are now thin glue +// over these; here we drive the logic on tiny problems with cheap analytic +// noise, and (with HDF5) round-trip the writers to pin the dataset-key contract +// with the Python analysis scripts. + +using namespace qew; + +namespace { + +struct ZeroNoise { + KOKKOS_INLINE_FUNCTION + NoiseSample sample(real_t, real_t) const { return NoiseSample{0, 0, 0}; } +}; + +// Washboard V = A cos(k y), no x dependence: a flat line is a single particle +// with analytic depinning threshold f_c = A k. +struct SinusoidNoise { + real_t A = 0.1; + real_t k = 2.0 * 3.14159265358979323846; // lambda = 1 + KOKKOS_INLINE_FUNCTION + NoiseSample sample(real_t, real_t y) const { + return NoiseSample{A * Kokkos::cos(k * y), 0, -A * k * Kokkos::sin(k * y)}; + } +}; + +} // namespace + +// ---- fmt(): the HDF5 key formatter -------------------------------------- +TEST(Fmt, MatchesPrintfSpecs) { + EXPECT_EQ(fmt("%.1g", 0.0), "0"); + EXPECT_EQ(fmt("%.3g", 0.1), "0.1"); + EXPECT_EQ(fmt("%.3g", 100.0), "100"); + EXPECT_EQ(fmt("%.6g", 1.0), "1"); + EXPECT_EQ(fmt("%.4g", 0.8), "0.8"); +} + +// ---- static_sweep ------------------------------------------------------- +TEST(Executables, StaticSweepShapeAndMapping) { + StaticConfig cfg; + cfg.model = Model::Linear; + cfg.amplitude = 1.3; + cfg.xi = 0.1; + cfg.Lx = 1.0; + cfg.Ly = 1.0; + cfg.nx = 16; + cfg.ny = 8; + cfg.pinning_lengths = {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()); + 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_NEAR(prof.line_tension, + std::pow(cfg.pinning_lengths[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); + const double span = *std::max_element(prof.h.begin(), prof.h.end()) - + *std::min_element(prof.h.begin(), prof.h.end()); + EXPECT_LT(span, 1e-9); + for (double v : prof.h) EXPECT_NEAR(v, 0.5 * cfg.Ly, 1e-9); + } +} + +// ---- dynamic_anneal ----------------------------------------------------- +TEST(Executables, DynamicAnnealTimestepAndSnapshots) { + DynamicConfig cfg; + cfg.model = Model::Linear; + cfg.amplitude = 1.0; + cfg.xi = 0.1; + cfg.pinning_length = 0.01; + cfg.Lx = 1.0; + cfg.Ly = 1.0; + cfg.nx = 32; + cfg.ny = 8; + cfg.tau = 1.0; + cfg.tl_start = 0.01; + cfg.tl_end = 0.0001; + cfg.total_time = 200.0; + + const DynamicResult res = dynamic_anneal(cfg, ZeroNoise{}); + + // dt = min(0.1, tau/10 * min(fac, fac^2)), fac = xi / max(tl_start, tl_end). + const real_t fac = cfg.xi / std::max(cfg.tl_start, cfg.tl_end); + 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); + EXPECT_NEAR(res.line_tension, + std::pow(cfg.pinning_length / cfg.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); + EXPECT_NEAR(res.snapshots.front().t, 10 * res.dt, res.dt); + for (std::size_t i = 1; i < res.snapshots.size(); ++i) { + EXPECT_GT(res.snapshots[i].t, res.snapshots[i - 1].t); + EXPECT_NEAR(res.snapshots[i].t / res.snapshots[i - 1].t, 10.0, 0.5); + } + for (const auto &snap : res.snapshots) { + ASSERT_EQ(snap.h.size(), static_cast(cfg.nx)); + for (double v : snap.h) EXPECT_TRUE(std::isfinite(v)); + } + EXPECT_LE(res.snapshots.back().t, res.total_time); +} + +// ---- depinning_ramp ----------------------------------------------------- +TEST(Executables, DepinningRampBracketsCriticalForce) { + const SinusoidNoise noise; + const real_t fc_exact = noise.A * noise.k; // = 0.2*pi ~ 0.628 + const real_t lambda = 2 * 3.14159265358979323846 / noise.k; + + DepinningConfig cfg; + cfg.model = Model::Linear; + cfg.amplitude = 1.0; + cfg.pinning_length = 1.0; + cfg.Lx = 8.0; + cfg.Ly = 1.0; + cfg.nx = 8; + cfg.ny = 8; + cfg.f_step = 0.2; + cfg.f_max = 2.0; + cfg.rp.dstep = lambda / 100; + cfg.rp.max_advance = lambda; + cfg.rp.runaway = 3 * lambda; + cfg.rp.root_tol = 1e-11; + cfg.rp.sweep_tol = 1e-8; + cfg.rp.max_sweeps = 20000; + + const DepinningResult res = depinning_ramp(cfg, noise); + + ASSERT_FALSE(res.steps.empty()); + EXPECT_TRUE(res.steps.front().result.blocked); // f = 0 blocks + EXPECT_EQ(res.steps.front().f, 0.0); + + // The ramp stops at the first unblocked force, which brackets f_c from above. + const DepinningStep &last = res.steps.back(); + EXPECT_FALSE(last.result.blocked); + EXPECT_TRUE(last.result.runaway); + EXPECT_TRUE(last.h.empty()); // no profile stored for the runaway force + EXPECT_NEAR(res.f_c_upper_bound, last.f, 1e-9); + + EXPECT_GT(res.f_c_upper_bound, 0.0); + EXPECT_LE(res.f_c_upper_bound, cfg.f_max); + // f_c lies in (f_c_upper_bound - f_step, f_c_upper_bound]. + EXPECT_LE(res.f_c_upper_bound - cfg.f_step, fc_exact); + EXPECT_LE(fc_exact, res.f_c_upper_bound + 1e-9); + + // Every stored (blocked) step carries a full profile; the runaway one alone + // is profile-free. + for (const auto &st : res.steps) { + if (st.result.blocked) + EXPECT_EQ(st.h.size(), static_cast(cfg.nx)); + } +} + +#ifdef QEW_WITH_HDF5 +namespace { + +bool has_dataset(HighFive::File &f, const std::string &path) { + return f.exist(path) && f.getObjectType(path) == HighFive::ObjectType::Dataset; +} + +} // namespace + +// The writers' dataset KEYS are a contract with the Python analysis scripts. +// Round-trip each writer and assert the exact key strings + values. +TEST(Executables, WriteStaticKeysAndValues) { + StaticConfig cfg; + cfg.Lx = 4.0; + cfg.driving_force = 0.0; + + std::vector profiles(1); + profiles[0].pinning_length = 0.1; + profiles[0].line_tension = 0.5; + profiles[0].h = {1.0, 2.0, 3.0, 4.0}; + + const std::string fname = "test_write_static.h5"; + { + HighFive::File file(fname, HighFive::File::Truncate); + write_static(file, cfg, profiles); + } + HighFive::File file(fname, HighFive::File::ReadOnly); + + double Lx = 0; + file.getDataSet("physical_size").read(Lx); + EXPECT_DOUBLE_EQ(Lx, cfg.Lx); + + const std::string prefix = "driving_force=0/pinning_length=0.1"; + ASSERT_TRUE(has_dataset(file, prefix + "/h")); + ASSERT_TRUE(has_dataset(file, prefix + "/pinning_length")); + std::vector h; + file.getDataSet(prefix + "/h").read(h); + EXPECT_EQ(h, profiles[0].h); + double lp = 0; + file.getDataSet(prefix + "/pinning_length").read(lp); + EXPECT_DOUBLE_EQ(lp, 0.1); +} + +TEST(Executables, WriteDynamicKeysAndValues) { + DynamicConfig cfg; + cfg.Lx = 2.0; + cfg.driving_force = 0.0; + cfg.pinning_length = 0.01; + + DynamicResult res; + res.line_tension = 2.0; + res.snapshots.push_back(DynamicSnapshot{static_cast(1.0), {5.0, 6.0}}); + + const std::string fname = "test_write_dynamic.h5"; + { + HighFive::File file(fname, HighFive::File::Truncate); + write_dynamic(file, cfg, res); + } + HighFive::File file(fname, HighFive::File::ReadOnly); + + const std::string prefix = "timestep=1/driving_force=0/line_tension=2"; + ASSERT_TRUE(has_dataset(file, prefix + "/h")); + ASSERT_TRUE(has_dataset(file, prefix + "/pinning_length")); + std::vector h; + file.getDataSet(prefix + "/h").read(h); + EXPECT_EQ(h, (std::vector{5.0, 6.0})); +} + +TEST(Executables, WriteDepinningSkipsRunawayAndWritesFc) { + DepinningConfig cfg; + cfg.Lx = 8.0; + cfg.pinning_length = 5.0; + + DepinningResult res; + res.f_c_upper_bound = 0.8; + DepinningStep blocked; + blocked.f = 0.0; + blocked.result.blocked = true; + blocked.h = {0.5, 0.5, 0.5}; + DepinningStep runaway; + runaway.f = 0.8; + runaway.result.blocked = false; + runaway.result.runaway = true; // h intentionally empty + res.steps = {blocked, runaway}; + + const std::string fname = "test_write_depinning.h5"; + { + HighFive::File file(fname, HighFive::File::Truncate); + write_depinning(file, cfg, res); + } + HighFive::File file(fname, HighFive::File::ReadOnly); + + ASSERT_TRUE(has_dataset(file, "driving_force=0/h")); + ASSERT_TRUE(has_dataset(file, "driving_force=0/pinning_length")); + // The runaway force stores no profile. + EXPECT_FALSE(file.exist("driving_force=0.8/h")); + + double fc = 0; + ASSERT_TRUE(has_dataset(file, "f_c_upper_bound")); + file.getDataSet("f_c_upper_bound").read(fc); + EXPECT_DOUBLE_EQ(fc, 0.8); +} +#endif // QEW_WITH_HDF5 diff --git a/tests/test_filtered_noise.cpp b/tests/test_filtered_noise.cpp index b625982..9c89006 100644 --- a/tests/test_filtered_noise.cpp +++ b/tests/test_filtered_noise.cpp @@ -229,6 +229,57 @@ TEST(FilteredNoise, GradientConsistentWithModel) { } } +// The pinning field is periodic in BOTH directions, and sample() must wrap any +// query coordinate -- including negative ones and many periods out -- back into +// the fundamental cell. The depinning and dynamics runs drive h forward through +// many periods of Ly (and the noise x is the edge-midpoint, which the periodic +// h advances past Lx), so correct wrapping is load-bearing for the physics. No +// other test queries outside [0,Lx) x [0,Ly); this pins it directly. +TEST(FilteredNoise, PeriodicWrapAcrossManyCells) { + const int nx = 48, ny = 40; + const real_t Lx = 2.5, Ly = 1.3; + FilteredNoise noise(nx, ny, Lx, Ly, 1.0, 0.25, 0.25, /*seed=*/17); + const DeviceNoise dn = noise.device_noise(); + + // Base query points strictly inside the cell. + std::vector xs, ys; + std::mt19937 rng(123); + std::uniform_real_distribution ux(0, Lx), uy(0, Ly); + const int m = 60; + for (int k = 0; k < m; ++k) { + xs.push_back(ux(rng)); + ys.push_back(uy(rng)); + } + std::vector v0, dx0, dy0; + sample_points(dn, xs, ys, v0, dx0, dy0); + + // Same points shifted by integer numbers of periods (including negative and + // large offsets) must reproduce value AND both derivatives bit-for-tight. + const std::pair shifts[] = { + {1, 0}, {0, 1}, {-1, 0}, {0, -1}, {3, -2}, {-5, 7}, {-13, -11}}; + for (const auto &sh : shifts) { + std::vector xss = xs, yss = ys; + for (int k = 0; k < m; ++k) { + xss[k] += sh.first * Lx; + yss[k] += sh.second * Ly; + } + std::vector v, dvx, dvy; + sample_points(dn, xss, yss, v, dvx, dvy); + double mv = 0, mdx = 0, mdy = 0; + for (int k = 0; k < m; ++k) { + mv = std::max(mv, std::abs(v[k] - v0[k])); + mdx = std::max(mdx, std::abs(dvx[k] - dx0[k])); + mdy = std::max(mdy, std::abs(dvy[k] - dy0[k])); + } + EXPECT_LT(mv, 1e-9) << "value not periodic for shift (" << sh.first + << "," << sh.second << ")"; + EXPECT_LT(mdx, 1e-9) << "dv/dx not periodic for shift (" << sh.first + << "," << sh.second << ")"; + EXPECT_LT(mdy, 1e-9) << "dv/dy not periodic for shift (" << sh.first + << "," << sh.second << ")"; + } +} + #ifdef QEW_WITH_HDF5 // Test 3c (golden): prefilter Python's noise_on_grid and check the interpolated // value/derivatives against SciPy's cubic B-spline at the same query points. diff --git a/tests/test_fire.cpp b/tests/test_fire.cpp index 953ca08..2952453 100644 --- a/tests/test_fire.cpp +++ b/tests/test_fire.cpp @@ -159,3 +159,30 @@ TEST(Fire, Deterministic) { const auto a = to_host(h1), b = to_host(h2); for (int i = 0; i < n; ++i) EXPECT_EQ(a[i], b[i]); } + +// Capped iterations on a problem that cannot converge in time must report +// converged=false with iterations==max_iter and a non-zero residual force -- +// the failure path callers branch on. Every other FIRE test asserts the success +// path, so a regression that always returned converged=true would slip through. +TEST(Fire, ReportsNonConvergence) { + const int n = 64; + Params p; + p.physical_size = n; // dx = 1 + p.line_tension = 1.0; + p.driving_force = 0.0; + p.model = Model::Linear; + HarmonicWell well{2.0, 0.5}; + + std::vector init(n, 0.0); // far from the minimum + View1D h = make_line(init); + + FireParams fp; + fp.dt_max = 0.1 * std::sqrt(1.0 / p.line_tension); + fp.ftol = 1e-12; // unreachable in a couple of steps + fp.max_iter = 3; + const FireResult r = fire_minimize(h, p, well, fp); + + EXPECT_FALSE(r.converged); + EXPECT_EQ(r.iterations, fp.max_iter); + EXPECT_GE(r.max_force, fp.ftol); +} diff --git a/tests/test_lbfgs.cpp b/tests/test_lbfgs.cpp index fee7676..440e6a0 100644 --- a/tests/test_lbfgs.cpp +++ b/tests/test_lbfgs.cpp @@ -101,6 +101,132 @@ TEST(Lbfgs, ConvergesToHarmonicMinimum) { EXPECT_LT(max_dev, 1e-6); } +// L-BFGS on the ARCLENGTH model: a flat line in a harmonic trough at y0 is still +// the unique minimiser (the arclength elastic term vanishes for a flat line), so +// it converges to h = y0. Every other L-BFGS test uses the linear model; this +// exercises the arclength branch of objective()/gradient() through the solver. +TEST(Lbfgs, ConvergesToHarmonicMinimumArclength) { + const int n = 96; + Params p = linear_params(n); + p.model = Model::Arclength; + HarmonicWell well{2.0, 0.5}; + View1D h = make_line(std::vector(n, 0.0)); + + LbfgsParams lp; + lp.ftol = 1e-9; + const LbfgsResult r = lbfgs_minimize(h, p, well, lp); + + EXPECT_TRUE(r.converged); + const auto hf = to_host(h); + real_t max_dev = 0; + for (real_t v : hf) max_dev = std::max(max_dev, std::abs(v - well.y0)); + EXPECT_LT(max_dev, 1e-6); +} + +// Capped iterations report converged=false (the failure path callers branch on). +// All other L-BFGS tests assert the success path. +TEST(Lbfgs, ReportsNonConvergence) { + const int n = 96; + Params p = linear_params(n); + HarmonicWell well{2.0, 0.5}; + View1D h = make_line(std::vector(n, 0.0)); + + LbfgsParams lp; + lp.ftol = 1e-12; // unreachable in one iteration + lp.max_iter = 1; + const LbfgsResult r = lbfgs_minimize(h, p, well, lp); + + EXPECT_FALSE(r.converged); + EXPECT_GE(r.max_force, lp.ftol); +} + +// Curvature safeguard (the `sy <= eps |s||y|` skip branch). By Cauchy-Schwarz +// sy <= |s||y| = sqrt(ss*yy), so an enormous curvature_eps makes the acceptance +// test sy > eps*sqrt(ss*yy) fail for EVERY pair: no (s,y) is ever stored. The +// solver must then degrade gracefully to (scaled) steepest descent and still +// reach the minimum -- the safeguard keeps the implicit inverse Hessian positive +// definite rather than corrupting the solve. A convex problem never skips on its +// own, so this is the only way to exercise the branch. +TEST(Lbfgs, CurvatureSkipForcedByLargeEps) { + const int n = 16; + Params p = linear_params(n); // small, well-conditioned -> SD converges + HarmonicWell well{2.0, 0.5}; + View1D h = make_line(std::vector(n, 0.0)); + + LbfgsParams lp; + lp.ftol = 1e-8; + lp.max_iter = 20000; + lp.curvature_eps = 1e30; // reject every curvature pair + const LbfgsResult r = lbfgs_minimize(h, p, well, lp); + + EXPECT_GT(r.n_curvature_skips, 0); + EXPECT_TRUE(r.converged); + const auto hf = to_host(h); + real_t max_dev = 0; + for (real_t v : hf) max_dev = std::max(max_dev, std::abs(v - well.y0)); + EXPECT_LT(max_dev, 1e-6); +} + +// Line-search give-up branch: max_ls = 0 means the backtracking loop runs zero +// trials, so the very first (steepest-descent) step "fails". With no history to +// fall back on, the solver reports non-convergence immediately (iteration 0) +// rather than looping or dereferencing empty history. +TEST(Lbfgs, LineSearchGivesUpWithNoHistory) { + const int n = 32; + Params p = linear_params(n); + HarmonicWell well{2.0, 0.5}; + View1D h = make_line(std::vector(n, 0.0)); + + LbfgsParams lp; + lp.ftol = 1e-9; + lp.max_ls = 0; + const LbfgsResult r = lbfgs_minimize(h, p, well, lp); + + EXPECT_FALSE(r.converged); + EXPECT_EQ(r.iterations, 0); + EXPECT_EQ(r.n_ls_restarts, 0); // never got far enough to build history +} + +// Line-search RESTART branch: when a unit L-BFGS step fails the line search but +// history exists, the solver discards its history and retries as steepest +// descent instead of giving up or dereferencing empty history. +// +// Trigger: a rough, strong-noise landscape on which a unit L-BFGS step +// overshoots a noise minimum by ~16x -- it needs ~5 Armijo backtracks to +// recover -- combined with max_ls clamped to 3. Once history is built, that +// step exhausts its 3 trials and falls into the restart branch. The overshoot +// is decisive (mls=1,2,3 all fire it; the line search only stops failing at +// mls>=5), so the branch is reached with margin rather than on an FP knife-edge. +TEST(Lbfgs, LineSearchRestartsOnStepFailure) { + const int nx = 64, ny = 64; + const real_t Lx = 3.2, Ly = 1.0; + FilteredNoise noise(nx, ny, Lx, Ly, /*amplitude=*/2.0, 0.15, 0.15, /*seed=*/3); + + Params p; + p.physical_size = Lx; + p.line_tension = 0.5; + p.driving_force = 0.0; + p.model = Model::Linear; + + View1D h = make_line(std::vector(nx, 0.5 * Ly)); + LbfgsParams lp; + lp.ftol = 1e-6; + lp.max_iter = 50000; + lp.max_ls = 3; // smaller than the ~5 backtracks the overshoot needs + const LbfgsResult r = lbfgs_minimize(h, p, noise.device_noise(), lp); + EXPECT_GT(r.n_ls_restarts, 0); // the restart branch executed + + // Control: with the normal backtracking budget the SAME problem never needs + // a restart and converges -- confirming the restart is induced by the clamp + // (a genuine safety net) and not a routine code path. + View1D h2 = make_line(std::vector(nx, 0.5 * Ly)); + LbfgsParams lp2 = lp; + lp2.max_ls = 30; + const LbfgsResult r2 = lbfgs_minimize(h2, p, noise.device_noise(), lp2); + EXPECT_TRUE(r2.converged); + EXPECT_EQ(r2.n_ls_restarts, 0); +} + // On a quadratic problem L-BFGS converges in few iterations (vs FIRE's many). TEST(Lbfgs, FastOnQuadratic) { const int n = 128; diff --git a/tests/test_qew_model.cpp b/tests/test_qew_model.cpp index 0c7513c..8cf99ca 100644 --- a/tests/test_qew_model.cpp +++ b/tests/test_qew_model.cpp @@ -33,6 +33,14 @@ struct AnalyticNoise { } }; +// Spatially constant pinning potential V == v0 (zero derivatives). Lets us write +// the total energy of a flat line in closed form, independent of the assembly. +struct ConstNoise { + real_t v0 = 0; + KOKKOS_INLINE_FUNCTION + NoiseSample sample(real_t, real_t) const { return NoiseSample{v0, 0, 0}; } +}; + void check_gradient(Model model) { const int n = 64; @@ -90,3 +98,66 @@ TEST(QewModel, GradientMatchesFiniteDifferenceLinear) { TEST(QewModel, GradientMatchesFiniteDifferenceArclength) { check_gradient(Model::Arclength); } + +// Absolute energy of a FLAT line in a constant potential is known in closed +// form: every edge slope is zero (no line-tension energy in either model), so +// E = sum_i dx * (v0 - f h_i) = L * (v0 - f c). +// This pins the absolute scale of objective() -- the dx weighting, the inclusion +// of the noise value, and the SIGN of the driving-force term -- for BOTH models +// (the gradient-consistency test above is blind to all three). The linear case +// is otherwise only checked indirectly via the SciPy L-BFGS golden; the +// arclength energy value is checked nowhere else. +TEST(QewModel, ObjectiveFlatLineConstantPotential) { + const int n = 50; + const real_t L = 3.0; + const real_t c = 0.7; // flat line height + const real_t v0 = 1.25; // constant potential value + const real_t f = 0.4; // driving force + + View1D h("h", n); + Kokkos::deep_copy(h, c); + + for (const Model model : {Model::Linear, Model::Arclength}) { + Params p; + p.physical_size = L; + p.line_tension = 2.3; // irrelevant for a flat line, but should be ignored + p.driving_force = f; + p.model = model; + + const real_t e = objective(h, p, ConstNoise{v0}); + const real_t expected = L * (v0 - f * c); + EXPECT_NEAR(e, expected, 1e-12) + << "model=" << (model == Model::Linear ? "linear" : "arclength"); + } +} + +// A single-mode line h_i = a sin(2 pi i / n) in the LINEAR model with zero noise +// and zero force has a closed-form discrete elastic energy: +// E = sum_i dx * (1/2) lt ((h_{i+1}-h_i)/dx)^2 +// = (lt a^2 n / (2 dx)) * (1 - cos(2 pi / n)) [periodic, exact]. +// Pins the linear line-tension term's value (not just its gradient) on a +// non-trivial profile. +TEST(QewModel, ObjectiveSinusoidLinearElastic) { + const int n = 32; + const real_t L = static_cast(n); // dx = 1 + const real_t dx = L / n; + const real_t lt = 1.7; + const real_t a = 0.9; + const real_t two_pi = 2 * 3.14159265358979323846; + + View1D h("h", n); + auto hh = Kokkos::create_mirror_view(h); + for (int i = 0; i < n; ++i) hh(i) = a * std::sin(two_pi * i / n); + Kokkos::deep_copy(h, hh); + + Params p; + p.physical_size = L; + p.line_tension = lt; + p.driving_force = 0.0; + p.model = Model::Linear; + + const real_t e = objective(h, p, ConstNoise{0.0}); + const real_t expected = + (lt * a * a * n / (2 * dx)) * (1 - std::cos(two_pi / n)); + EXPECT_NEAR(e, expected, 1e-10); +} diff --git a/tests/test_rosso_krauth.cpp b/tests/test_rosso_krauth.cpp index 9eeae07..8d2c828 100644 --- a/tests/test_rosso_krauth.cpp +++ b/tests/test_rosso_krauth.cpp @@ -213,3 +213,78 @@ TEST(RossoKrauth, Deterministic) { const auto a = to_host(h1), b = to_host(h2); for (int i = 0; i < nx; ++i) EXPECT_EQ(a[i], b[i]); } + +// Above threshold the construction must explicitly report RUNAWAY (not merely +// "not blocked"): the line advances past rp.runaway with no metastable state. +// The analytic-fc bisection test only ever reads `.blocked`, so the `runaway` +// flag and the runaway exit branch are otherwise unexercised -- yet the +// depinning executable keys its f_c bracketing and its "DEPINNED" report on it. +TEST(RossoKrauth, RunawayFlagSetAboveThreshold) { + const int n = 16; + const SinusoidNoise noise; + const real_t fc_exact = noise.A * noise.k; + const real_t lambda = 2 * 3.14159265358979323846 / noise.k; + + Params p; + p.physical_size = n; // dx = 1 (flat line, elastic irrelevant) + p.line_tension = 1.0; + p.driving_force = 2.0 * fc_exact; // well above threshold -> runs away + p.model = Model::Linear; + + RkParams rp; + rp.dstep = lambda / 200; + rp.max_advance = lambda; + rp.runaway = 3 * lambda; + rp.root_tol = 1e-12; + rp.sweep_tol = 1e-11; + rp.max_sweeps = 100000; + + View1D h = flat_line(n, 0.0); + const RkResult r = rk_block(h, p, noise, rp); + EXPECT_TRUE(r.runaway); + EXPECT_FALSE(r.blocked); + EXPECT_GT(r.mean_h, rp.runaway); // advanced past the runaway threshold +} + +// Below threshold the arclength model also blocks into an arrested state. The +// analytic-fc and forward-and-arrested tests pin only the LINEAR model through +// rk_block (Arclength appears only in the determinism test, which checks +// reproducibility, not correctness). Here the arclength site velocity must be +// <= 0 everywhere at the returned configuration. +TEST(RossoKrauth, BlocksArclengthModel) { + const int nx = 64, ny = 64; + const real_t Lx = 3.2, Ly = 1.0; + FilteredNoise noise(nx, ny, Lx, Ly, /*amplitude=*/0.6, 0.2, 0.2, /*seed=*/29); + const DeviceNoise dn = noise.device_noise(); + + Params p; + p.physical_size = Lx; + p.line_tension = 1.2; + p.driving_force = 0.1; + p.model = Model::Arclength; + + View1D h = flat_line(nx, 0.5 * Ly); + + RkParams rp; + rp.dstep = 0.5 * (Ly / ny); + rp.max_advance = Ly; + rp.runaway = 3 * Ly; + rp.root_tol = 1e-11; + rp.sweep_tol = 1e-9; + rp.max_sweeps = 200000; + + const RkResult r = rk_block(h, p, dn, rp); + EXPECT_TRUE(r.blocked); + + View1D grad("grad", nx); + gradient(h, p, dn, grad); + real_t max_v = 0; + Kokkos::parallel_reduce( + "max_v", nx, + KOKKOS_LAMBDA(const int i, real_t &acc) { + const real_t v = -grad(i); + if (v > acc) acc = v; + }, + Kokkos::Max(max_v)); + EXPECT_LT(max_v, 1e-6); +}