diff --git a/compute.py b/compute.py index aef23cc..ddb5d6b 100644 --- a/compute.py +++ b/compute.py @@ -894,6 +894,7 @@ def run_engine(engine, input_dir, output_dir, config): "damping_force": int(config.get("damping_force", 0)), "gravity_strength": float(config.get("gravity_strength", 1.0)), "driving_force": int(config.get("driving_force", 0)), + "save_trajectory": int(config.get("save_trajectory", 0)), } param_path = os.path.join(script_dir, "engines", engine, "param.json") os.makedirs(os.path.dirname(param_path), exist_ok=True) @@ -1008,7 +1009,17 @@ def run_engine(engine, input_dir, output_dir, config): except OSError: pass - # 加载输出的 trajectory.txt + # 加载输出的 trajectory.txt / display.txt + save_traj = int(config.get("save_trajectory", 0)) + if not save_traj: + # save_trajectory=0：引擎只写 display.txt + disp_path = os.path.join(os.path.abspath(output_dir), "display.txt") + if not os.path.exists(disp_path): + raise FileNotFoundError(f"引擎未生成 display.txt: {disp_path}") + print(f"[compute] 引擎已生成 {disp_path}") + return None, None, None, None, None, None + + # save_trajectory=1：加载完整 trajectory.txt traj_path = os.path.join(os.path.abspath(output_dir), "trajectory.txt") if not os.path.exists(traj_path): raise FileNotFoundError(f"引擎未生成 trajectory.txt: {traj_path}") diff --git a/dynamics.py b/dynamics.py index a6d59f3..d601b51 100644 --- a/dynamics.py +++ b/dynamics.py @@ -246,15 +246,15 @@ def run_case(config_path, runtime_base, input_dir="input", output_dir="output", else: print("[run] 步骤 [模拟] 已跳过") - # 3. 生成 display.txt：模拟完成后总是刷新，避免参数变更后缓存过时 + # 3. 检查/生成 display.txt disp_path_new = os.path.join(output_dir_abs, "display.txt") + save_traj = int(config.get("save_trajectory", 0)) - if engine == "python": - # Python 引擎在 run_simulation 内部已写入 display.txt - if os.path.exists(disp_path_new): - print(f"[run] 发现已有 display.txt（引擎直接抽帧）") - else: - # 外部引擎：从 trajectory.txt 重新抽帧（覆盖旧的 display.txt） + if os.path.exists(disp_path_new): + # Python 引擎或新版外部引擎（save_trajectory=0）已直接写入 + print(f"[run] 发现已有 display.txt（引擎直接抽帧）") + elif engine != "python" and os.path.exists(os.path.join(output_dir_abs, "trajectory.txt")): + # 旧版外部引擎：从 trajectory.txt 抽帧 traj_path = os.path.join(output_dir_abs, "trajectory.txt") if not os.path.exists(traj_path): print(f"[run] 错误: 找不到 trajectory.txt 或 display.txt") @@ -310,15 +310,13 @@ def run_case(config_path, runtime_base, input_dir="input", output_dir="output", header_fields=hf) print(f"[run] 从 trajectory.txt 抽帧生成 display.txt ({n_frames} 帧)") - # 外部引擎：save_trajectory=0 时清理 trajectory.txt - save_traj = int(config.get("save_trajectory", 0)) - if engine != "python" and not save_traj: - traj_path = os.path.join(output_dir_abs, "trajectory.txt") - try: - os.remove(traj_path) - print(f"[run] save_trajectory=0，已删除 {traj_path}") - except OSError: - pass + # save_trajectory=0 时清理 trajectory.txt + if not save_traj: + try: + os.remove(traj_path) + print(f"[run] save_trajectory=0，已删除 {traj_path}") + except OSError: + pass # 4. 绘图（可选） if not no_plot and config.get("step_plot", 1): diff --git a/engines/c/main.c b/engines/c/main.c index 61f4468..4e4fc7c 100644 --- a/engines/c/main.c +++ b/engines/c/main.c @@ -38,6 +38,7 @@ typedef struct { int damping_force; /* 阻尼开关 */ double gravity_strength; /* 万有引力强度 */ int driving_force; /* 驱动力开关 */ + int save_trajectory; /* 是否保存完整轨迹文件 */ } SimParams; /* ======================================================================== @@ -304,6 +305,7 @@ static SimParams read_params(const char *path) { p.damping_force = json_read_int(buf, "damping_force"); p.gravity_strength = json_read_double(buf, "gravity_strength"); p.driving_force = json_read_int(buf, "driving_force"); + p.save_trajectory = json_read_int(buf, "save_trajectory"); g_gravity_field = p.gravity_field; g_gravity_interaction = p.gravity_interaction; g_elastic_force = p.elastic_force; @@ -858,6 +860,49 @@ static void write_trajectory_json(const char *path, const Trajectory *traj, fclose(f); } +static void write_display_txt(const char *path, const Trajectory *traj, + const SimParams *params, const AtomData *atoms) +{ + FILE *f = fopen(path, "w"); + if (!f) die("无法写入 display.txt"); + + int n_frames = traj->n_steps; + int n_particles = traj->n_atoms; + int dynamic_steps = params->NT - params->warmup_steps; + double T_total = dynamic_steps * params->DT; + + fprintf(f, "number of frames: %d\n", n_frames); + fprintf(f, "number of particles: %d\n", n_particles); + fprintf(f, "DT: %.16g\n", params->DT); + fprintf(f, "NSTEP: %d\n", params->NSTEP); + fprintf(f, "method: %s\n", params->method); + fprintf(f, "warmup_steps: %d\n", params->warmup_steps); + fprintf(f, "dynamic_steps: %d\n", dynamic_steps); + fprintf(f, "T_total: %.16g\n", T_total); + fprintf(f, "box_a: %.16g\n", params->box_a); + + if (params->driving_force) { + fprintf(f, "driving_force: 1\n"); + } else { + fprintf(f, "driving_force: 0\n"); + } + fprintf(f, "\n"); + + for (int t = 0; t < n_frames; t++) { + fprintf(f, "frame: %3d\n", t + 1); + fprintf(f, "n x y z vx vy vz\n"); + for (int i = 0; i < n_particles; i++) { + int idx = t * n_particles + i; + fprintf(f, "%4d %12.6f %12.6f %12.6f %10.6f %10.6f %10.6f\n", + atoms->atom_ids[i], + traj->x[idx], traj->y[idx], traj->z[idx], + traj->vx[idx], traj->vy[idx], traj->vz[idx]); + } + } + + fclose(f); +} + // ======================================================================== // 主函数 // ======================================================================== @@ -902,17 +947,29 @@ int main(int argc, char **argv) { vz[i] = atoms.vel_0[i*3+2]; } - /* 分配轨迹缓冲区：改用 record_steps */ + /* 分配轨迹缓冲区 */ int record_steps = params.NT - params.warmup_steps; Trajectory traj; - traj.n_steps = record_steps; traj.n_atoms = n; - traj.x = (double*)xmalloc(record_steps * n * sizeof(double) * 6); - traj.y = traj.x + record_steps * n; - traj.z = traj.y + record_steps * n; - traj.vx = traj.z + record_steps * n; - traj.vy = traj.vx + record_steps * n; - traj.vz = traj.vy + record_steps * n; + if (params.save_trajectory) { + traj.n_steps = record_steps; + traj.x = (double*)xmalloc(record_steps * n * sizeof(double) * 6); + traj.y = traj.x + record_steps * n; + traj.z = traj.y + record_steps * n; + traj.vx = traj.z + record_steps * n; + traj.vy = traj.vx + record_steps * n; + traj.vz = traj.vy + record_steps * n; + } else { + int sampled_steps = (record_steps + params.NSTEP - 1) / params.NSTEP; + if (sampled_steps < 1) sampled_steps = 1; + traj.n_steps = sampled_steps; + traj.x = (double*)xmalloc(sampled_steps * n * sizeof(double) * 6); + traj.y = traj.x + sampled_steps * n; + traj.z = traj.y + sampled_steps * n; + traj.vx = traj.z + sampled_steps * n; + traj.vy = traj.vx + sampled_steps * n; + traj.vz = traj.vy + sampled_steps * n; + } /* 预热 */ /* 初始时刻 t=0 驱动力（与 Python run_simulation 一致）*/ @@ -930,6 +987,7 @@ int main(int argc, char **argv) { /* 记录 */ int _prog_interval = record_steps / 100; if (_prog_interval < 1) _prog_interval = 1; + int sample_idx = 0; for (int s = 0; s < record_steps; s++) { if (s % _prog_interval == 0 && s > 0) { printf("[C-engine] progress: %d/%d\n", s, record_steps); @@ -937,13 +995,18 @@ int main(int argc, char **argv) { } double t = (s + params.warmup_steps) * params.DT; if (params.driving_force) apply_driving_force(n, x, y, z, vx, vy, vz, t, s, params.DT, &drivers); - for (int i = 0; i < n; i++) { - traj.x[ s * n + i] = x[i]; - traj.y[ s * n + i] = y[i]; - traj.z[ s * n + i] = z[i]; - traj.vx[s * n + i] = vx[i]; - traj.vy[s * n + i] = vy[i]; - traj.vz[s * n + i] = vz[i]; + int do_record = params.save_trajectory || (s % params.NSTEP == 0); + if (do_record) { + int idx = params.save_trajectory ? s : sample_idx; + for (int i = 0; i < n; i++) { + traj.x[ idx * n + i] = x[i]; + traj.y[ idx * n + i] = y[i]; + traj.z[ idx * n + i] = z[i]; + traj.vx[idx * n + i] = vx[i]; + traj.vy[idx * n + i] = vy[i]; + traj.vz[idx * n + i] = vz[i]; + } + if (!params.save_trajectory) sample_idx++; } apply_step(params.method, n, x, y, z, vx, vy, vz, atoms.masses, params.G, params.B, &bonds, atoms.fixed, @@ -952,8 +1015,13 @@ int main(int argc, char **argv) { } char out_path[512]; - snprintf(out_path, sizeof(out_path), "%s/trajectory.txt", output_dir); - write_trajectory_json(out_path, &traj, ¶ms, &atoms, &bonds); + if (params.save_trajectory) { + snprintf(out_path, sizeof(out_path), "%s/trajectory.txt", output_dir); + write_trajectory_json(out_path, &traj, ¶ms, &atoms, &bonds); + } else { + snprintf(out_path, sizeof(out_path), "%s/display.txt", output_dir); + write_display_txt(out_path, &traj, ¶ms, &atoms); + } clock_t t1 = clock(); double elapsed = (double)(t1 - t0) / CLOCKS_PER_SEC; diff --git a/engines/cpp/main.cpp b/engines/cpp/main.cpp index 6221cb2..0840213 100644 --- a/engines/cpp/main.cpp +++ b/engines/cpp/main.cpp @@ -43,6 +43,7 @@ struct SimParams { int damping_force = 0; double gravity_strength = 1.0; int driving_force = 0; + int save_trajectory = 1; }; // ======================================================================== @@ -131,6 +132,11 @@ static std::string json_read_string(const std::string &json, const std::string & return result; } +/* 检查 JSON 中是否存在某个 key */ +static bool json_has_key(const std::string &json, const std::string &key) { + return json.find("\"" + key + "\"") != std::string::npos; +} + /* 读取 JSON 数组 (如 "G": [0, 0, -9.8]) 到 double[3] */ static void json_read_double3(const std::string &json, const std::string &key, double out[3]) { auto pos = json.find("\"" + key + "\""); @@ -165,6 +171,9 @@ static SimParams read_params(const std::string &path) { p.damping_force = json_read_int(buf, "damping_force"); p.gravity_strength = json_read_double(buf, "gravity_strength"); p.driving_force = json_read_int(buf, "driving_force"); + // save_trajectory 默认 1（全量），仅在 JSON 中存在该 key 时覆盖 + if (json_has_key(buf, "save_trajectory")) + p.save_trajectory = json_read_int(buf, "save_trajectory"); return p; } @@ -681,6 +690,52 @@ static void apply_driving_force( } } +// ======================================================================== +// display.txt 输出（save_trajectory=0 时使用） +// ======================================================================== + +static void write_display_txt( + const std::string &path, + const std::vector &x, const std::vector &y, + const std::vector &z, const std::vector &vx, + const std::vector &vy, const std::vector &vz, + int n_steps, int n_atoms, + const SimParams ¶ms, const AtomData &atoms) +{ + std::ofstream f(path); + if (!f) die("无法写入 " + path); + std::cout << "[Cpp-engine] 正在写入显示数据…" << std::endl; + + double T_total = params.NT * params.DT; + + f << "number of frames: " << n_steps << "\n"; + f << "number of particles: " << n_atoms << "\n"; + f << "DT: " << params.DT << "\n"; + f << "NSTEP: " << params.NSTEP << "\n"; + f << "method: " << params.method << "\n"; + f << "warmup_steps: " << params.warmup_steps << "\n"; + f << "dynamic_steps: " << (params.NT - params.warmup_steps) << "\n"; + f << "T_total: " << T_total << "\n"; + f << "box_a: " << params.box_a << "\n"; + f << "driving_force: " << params.driving_force << "\n\n"; + + f << std::fixed << std::setprecision(6); + for (int t = 0; t < n_steps; t++) { + f << "frame: " << (t + 1) << "\n"; + f << "n x y z vx vy vz\n"; + for (int i = 0; i < n_atoms; i++) { + int base = t * n_atoms + i; + f << std::setw(3) << atoms.ids[i] + << std::setw(12) << x[base] + << std::setw(12) << y[base] + << std::setw(12) << z[base] + << std::setw(12) << vx[base] + << std::setw(12) << vy[base] + << std::setw(12) << vz[base] << "\n"; + } + } +} + // ======================================================================== // JSON 输出 // ======================================================================== @@ -817,12 +872,14 @@ int main(int argc, char **argv) { // 分配轨迹缓冲区 int record_steps = params.NT - params.warmup_steps; - std::vector traj_x(record_steps * n); - std::vector traj_y(record_steps * n); - std::vector traj_z(record_steps * n); - std::vector traj_vx(record_steps * n); - std::vector traj_vy(record_steps * n); - std::vector traj_vz(record_steps * n); + int nstep_sampling = (params.save_trajectory == 0) ? params.NSTEP : 1; + int buf_steps = (params.save_trajectory == 0) ? (record_steps / nstep_sampling) : record_steps; + std::vector traj_x(buf_steps * n); + std::vector traj_y(buf_steps * n); + std::vector traj_z(buf_steps * n); + std::vector traj_vx(buf_steps * n); + std::vector traj_vy(buf_steps * n); + std::vector traj_vz(buf_steps * n); // 预热 // 初始时刻 t=0 驱动力（与 Python run_simulation 一致） @@ -846,6 +903,7 @@ int main(int argc, char **argv) { // 记录 int _prog_int = record_steps / 100; if (_prog_int < 1) _prog_int = 1; + int si = 0; // 采样帧索引 for (int s = 0; s < record_steps; s++) { if (s % _prog_int == 0 && s > 0) { std::cout << "[Cpp-engine] progress: " << s << "/" << record_steps << std::endl; @@ -853,14 +911,17 @@ int main(int argc, char **argv) { double t = (s + params.warmup_steps) * params.DT; if (params.driving_force) apply_driving_force(n, x.data(), y.data(), z.data(), vx.data(), vy.data(), vz.data(), t, s, params.DT, drivers); - // 保存当前帧 - for (int i = 0; i < n; i++) { - traj_x[s * n + i] = x[i]; - traj_y[s * n + i] = y[i]; - traj_z[s * n + i] = z[i]; - traj_vx[s * n + i] = vx[i]; - traj_vy[s * n + i] = vy[i]; - traj_vz[s * n + i] = vz[i]; + // 保存当前帧（采样模式仅每 NSTEP 步保存一次） + if (s % nstep_sampling == 0) { + for (int i = 0; i < n; i++) { + traj_x[si * n + i] = x[i]; + traj_y[si * n + i] = y[i]; + traj_z[si * n + i] = z[i]; + traj_vx[si * n + i] = vx[i]; + traj_vy[si * n + i] = vy[i]; + traj_vz[si * n + i] = vz[i]; + } + si++; } apply_step(params.method, n, x.data(), y.data(), z.data(), @@ -874,9 +935,15 @@ int main(int argc, char **argv) { } // 输出轨迹 - std::string out_path = output_dir + "/trajectory.txt"; - write_trajectory_json(out_path, traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz, - record_steps, n, params, atoms, bonds); + if (params.save_trajectory == 0) { + std::string out_path = output_dir + "/display.txt"; + write_display_txt(out_path, traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz, + si, n, params, atoms); + } else { + std::string out_path = output_dir + "/trajectory.txt"; + write_trajectory_json(out_path, traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz, + record_steps, n, params, atoms, bonds); + } auto t1 = std::chrono::high_resolution_clock::now(); double elapsed = std::chrono::duration(t1 - t0).count();