Skip to content

test: cover periodic noise wrap, energy values, driving force, RK runaway, and minimizer failure paths #1

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
pastewka merged 3 commits into from
Jun 9, 2026
Merged

test: cover periodic noise wrap, energy values, driving force, RK runaway, and minimizer failure paths #1

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
4ce4ef5
test: cover periodic noise wrap, energy values, driving force, RK run…
claude Jun 8, 2026
cccfa57
refactor: extract executable orchestration + HDF5 I/O for testability
claude Jun 8, 2026
a453acd
test: cover L-BFGS curvature-skip and line-search restart/give-up bra…
claude Jun 9, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
118 changes: 40 additions & 78 deletions executables/qew_depinning/main.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 <cmath>
#include <cstdint>
#include <cstdio>
#include <string>
#include <vector>

#include <Kokkos_Core.hpp>
#include <highfive/H5File.hpp>

#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<real_t>(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<real_t>(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<real_t>(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<double>(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<double>(pinning_length),
static_cast<double>(line_tension));

View1D h("h", nx);
Kokkos::deep_copy(h, static_cast<real_t>(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<double>(cfg.pinning_length),
static_cast<double>(res.line_tension));
for (const DepinningStep &st : res.steps) {
std::printf("f=%7.3f sweeps=%-7d mean_h=%8.4f %s\n",
static_cast<double>(f), r.sweeps,
static_cast<double>(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<double> 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<double>(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<double>(f));
break;
}
static_cast<double>(st.f), st.result.sweeps,
static_cast<double>(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<double>(f_c));
if (res.f_c_upper_bound >= 0)
std::printf("Line depinned near f=%.3f; stopping ramp.\n",
static_cast<double>(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());
}
Expand Down
111 changes: 31 additions & 80 deletions executables/qew_dynamic/main.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 <algorithm>
#include <cmath>
#include <cstdint>
#include <cstdio>
#include <string>
#include <vector>

#include <Kokkos_Core.hpp>
#include <highfive/H5File.hpp>

#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<real_t>(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<real_t>(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<real_t>(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<double>(res.dt),
static_cast<double>(res.total_time),
static_cast<long>(res.total_time / res.dt),
static_cast<double>(res.line_tension));
for (const DynamicSnapshot &snap : res.snapshots)
std::printf(" saved t=%.6g\n", static_cast<double>(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<double>(Lx));

std::printf("dt=%.3g total_time=%.3g steps=%ld line_tension=%.4g\n",
static_cast<double>(dt), static_cast<double>(total_time),
static_cast<long>(total_time / dt),
static_cast<double>(line_tension));

auto save = [&](real_t t) {
auto hh = Kokkos::create_mirror_view(h);
Kokkos::deep_copy(hh, h);
std::vector<double> 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<double>(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());
Expand Down
Loading
Loading