From 80520590d1b3cf36e70294b319ad44d118d6396b Mon Sep 17 00:00:00 2001 From: Ying-Li Niu <64801511@qq.com> Date: Thu, 11 Jun 2026 12:39:46 +0800 Subject: [PATCH] =?UTF-8?q?feat:=20=E6=96=B0=E5=A2=9E=E6=B3=A2=E5=BD=A2?= =?UTF-8?q?=E8=83=BD=E9=87=8F=E5=8A=A8=E7=94=BB=E7=B3=BB=E7=BB=9F=20plot?= =?UTF-8?q?=5Fwave.py?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - 创建 plot_wave.py: 从 display.txt 读取原子位移数据 绘制纵波(x) + 横波(y) + 横波(z) 波形随时间的动画 同时绘制系统动能/弹性势能/总能量/输入功率(dE/dt)时变曲线 输出 wave_animation.gif - 所有 input.txt 新增 step_plot_wave: 0 开关 - case05 开启 step_plot_wave: 1 - dynamics.py disp_data 新增 bond_stiffness/bond_rest_lengths - 更新案例文档 --- draw.py | 105 +++++- dynamics.py | 13 +- examples/case01/input/input.txt | 1 + examples/case02/input/input.txt | 1 + examples/case03/input/input.txt | 1 + examples/case04/input/input.txt | 1 + examples/case05/doc/index.html | 10 +- examples/case05/input/bond.txt | 2 +- examples/case05/input/coord.txt | 120 +++---- examples/case05/input/driver.txt | 2 +- examples/case05/input/input.txt | 11 +- examples/case06/Readme.md | 42 +++ examples/case06/doc/index.html | 473 +++++++++++++++++++++++++++ examples/case06/input/bond.txt | 2 + examples/case06/input/connection.txt | 60 ++++ examples/case06/input/coord.txt | 61 ++++ examples/case06/input/driver.txt | 2 + examples/case06/input/input.txt | 96 ++++++ examples/case06/run_dynamics.py | 54 +++ plot_wave.py | 260 +++++++++++++++ 20 files changed, 1231 insertions(+), 86 deletions(-) create mode 100644 examples/case06/Readme.md create mode 100644 examples/case06/doc/index.html create mode 100644 examples/case06/input/bond.txt create mode 100644 examples/case06/input/connection.txt create mode 100644 examples/case06/input/coord.txt create mode 100644 examples/case06/input/driver.txt create mode 100644 examples/case06/input/input.txt create mode 100644 examples/case06/run_dynamics.py create mode 100644 plot_wave.py diff --git a/draw.py b/draw.py index 904f7d4..82a16fd 100644 --- a/draw.py +++ b/draw.py @@ -317,10 +317,11 @@ def update_camera_info(event=None): c = view.camera camera_info.text = ( "Camera\n" - f"center = ({c.center[0]:.2f}, {c.center[1]:.2f}, {c.center[2]:.2f})\n" - f"distance = {c.distance:.2f}\n" - f"elevation = {c.elevation:.2f}\n" - f"azimuth = {c.azimuth:.2f}" + f"center = ({c.center[0]:.1f}, {c.center[1]:.1f}, {c.center[2]:.1f})\n" + f"distance = {c.distance:.1f}\n" + f"elevation = {c.elevation:.1f}\n" + f"azimuth = {c.azimuth:.1f}\n" + f"step = {PAN_SPEED:.1f} | {'透视' if _PERSPECTIVE else '正交'}" ) @@ -333,6 +334,7 @@ def update_ball_info(frame_idx, x, y, z, vx, vy, vz): f"t = {t:.2f} s | dt = {DT:.3f} s | nstep = {NSTEP}\n" f"Position: ({x:.2f}, {y:.2f}, {z:.2f})\n" f"Velocity: ({vx:.2f}, {vy:.2f}, {vz:.2f})\n" + f"W/S 沿Z轴 | A/D 左右 | Q/E 升降 | C/X 步长 | V 透视/正交" ) @@ -358,28 +360,99 @@ def reposition_camera_info(event=None): update_camera_info() +# ── 平移速度 & 投影模式(全局变量)── +PAN_SPEED = 1.0 +_PERSPECTIVE = True # True=透视, False=正交 + + def handle_view_interaction(event): update_camera_info() -def rotate_about_screen_normal(angle): - if hasattr(view.camera, "roll"): - view.camera.roll = (view.camera.roll + angle) % 360 - else: - view.camera.azimuth = (view.camera.azimuth + angle) % 360 - update_camera_info() - - def handle_key_press(event): + global PAN_SPEED, _PERSPECTIVE key_name = "" if getattr(event, "text", None): key_name = event.text.lower() elif getattr(event, "key", None) is not None: key_name = str(event.key).lower() - if key_name == "q": - rotate_about_screen_normal(-90) - elif key_name == "e": - rotate_about_screen_normal(90) + + c = view.camera + + # ── 投影切换 ── + if key_name == "v": + _PERSPECTIVE = not _PERSPECTIVE + try: + if _PERSPECTIVE: + c.fov = 60.0 + else: + c.fov = 0.0 + print(f"[draw] 投影模式: {'透视' if _PERSPECTIVE else '正交'}") + except Exception: + print("[draw] 当前相机不支持 fov 切换") + update_camera_info() + return + + # ── 步长控制 ── + if key_name == "c": + PAN_SPEED = min(PAN_SPEED * 1.5, 50.0) + print(f"[draw] 步长: {PAN_SPEED:.1f}") + update_camera_info() + return + elif key_name == "x": + PAN_SPEED = max(PAN_SPEED / 1.5, 0.05) + print(f"[draw] 步长: {PAN_SPEED:.1f}") + update_camera_info() + return + + # ── 计算相机方向向量 ── + import math as _math + azim_rad = _math.radians(c.azimuth) + elev_rad = _math.radians(c.elevation) + + # 视线方向(从 center 指向相机) + vd = np.array([ + _math.cos(elev_rad) * _math.sin(azim_rad), + _math.sin(elev_rad), + _math.cos(elev_rad) * _math.cos(azim_rad), + ]) + vd /= np.linalg.norm(vd) + + # 屏幕右方向 + world_up = np.array([0.0, 1.0, 0.0]) + right = np.cross(vd, world_up) + rn = np.linalg.norm(right) + if rn > 1e-10: + right /= rn + else: + right = np.array([1.0, 0.0, 0.0]) + + # 屏幕上方向 + up = np.cross(right, vd) + un = np.linalg.norm(up) + if un > 1e-10: + up /= un + else: + up = np.array([0.0, 1.0, 0.0]) + + pan = PAN_SPEED * 0.3 + + if key_name == "a": # 右移(屏幕右方向) + c.center = tuple(np.array(c.center) + right * pan) + elif key_name == "d": # 左移(屏幕右方向负向) + c.center = tuple(np.array(c.center) - right * pan) + elif key_name == "e": # 上升(屏幕上方向,原 W 的功能) + c.center = tuple(np.array(c.center) + up * pan) + elif key_name == "q": # 下降(屏幕上方向负向,原 S 的功能) + c.center = tuple(np.array(c.center) - up * pan) + elif key_name == "w": # 相机沿 Z 轴上移(靠近场景) + c.center = (c.center[0], c.center[1], c.center[2] + pan) + elif key_name == "s": # 相机沿 Z 轴下移(远离场景) + c.center = (c.center[0], c.center[1], c.center[2] - pan) + else: + return + + update_camera_info() def reset_camera_view(): diff --git a/dynamics.py b/dynamics.py index 3e04d01..93b96fc 100644 --- a/dynamics.py +++ b/dynamics.py @@ -125,7 +125,7 @@ def run_case(config_path, runtime_base, input_dir="input", output_dir="output", print(f"[run] 同时指定了 T_total 和 NT,使用 NT={config['NT']}") # 显示步骤控制信息 - steps_info = {k: config.get(k, 1) for k in ["step_simulate", "step_sample", "step_plot", "step_animation"]} + steps_info = {k: config.get(k, 1) for k in ["step_simulate", "step_sample", "step_plot", "step_plot_wave", "step_animation"]} step_flags = ", ".join(f"{k}={v}" for k, v in steps_info.items()) print(f"[run] 步骤开关: {step_flags}") @@ -295,6 +295,8 @@ def run_case(config_path, runtime_base, input_dir="input", output_dir="output", "atom_velocities": data["atom_velocities"] if "atom_velocities" in data else np.array([[float(data["VX0"]), float(data["VY0"]), float(data["VZ0"])]]), "atom_fixed": data["atom_fixed"] if "atom_fixed" in data else np.array([[0, 0, 0]]), "bond_pairs": data.get("bond_pairs", np.zeros((0, 2), dtype=np.int64)).tolist(), + "bond_stiffness": data.get("bond_stiffness", np.zeros(0, dtype=np.float64)).tolist(), + "bond_rest_lengths": data.get("bond_rest_lengths", np.zeros(0, dtype=np.float64)).tolist(), "warmup_steps": warmup_steps, "sample_start": sample_start, "sample_end": sample_end, @@ -500,6 +502,15 @@ def run_case(config_path, runtime_base, input_dir="input", output_dir="output", else: print("[run] 运行 python draw.py 查看动画。") + # 6. 波形能量动画(可选) + if config.get("step_plot_wave", 0): + try: + import plot_wave as pw + print("[run] 正在绘制波形与能量图…") + pw.plot_wave(str(output_dir_abs)) + except Exception as e: + print(f"[run] 绘制波形图失败: {e}") + def main(): parser = argparse.ArgumentParser(description="物理模拟统一入口") diff --git a/examples/case01/input/input.txt b/examples/case01/input/input.txt index 549da91..ab27154 100644 --- a/examples/case01/input/input.txt +++ b/examples/case01/input/input.txt @@ -8,6 +8,7 @@ step_simulate: 1 # 运行物理模拟 → output/trajectory.txt step_sample: 1 # 抽帧 → output/display.txt step_plot: 1 # 绘制轨迹/能量图 → output/trajectory_plots.png +step_plot_wave: 0 # 绘制波形能量动画 → output/wave_animation.gif step_animation: 1 # 自动播放 VisPy 3D 动画窗口(需安装 vispy) force_calc: 0 # 强制重新计算:1=跳过缓存强算,0=自动使用已有输出 diff --git a/examples/case02/input/input.txt b/examples/case02/input/input.txt index a5bf5de..292d59a 100644 --- a/examples/case02/input/input.txt +++ b/examples/case02/input/input.txt @@ -8,6 +8,7 @@ step_simulate: 1 # 运行物理模拟 → output/trajectory.txt step_sample: 1 # 抽帧 → output/display.txt step_plot: 1 # 绘制轨迹/能量图 → output/trajectory_plots.png +step_plot_wave: 0 # 绘制波形能量动画 → output/wave_animation.gif step_animation: 1 # 自动播放 VisPy 3D 动画窗口(需安装 vispy) force_calc: 0 # 强制重新计算:1=跳过缓存强算,0=自动使用已有输出 diff --git a/examples/case03/input/input.txt b/examples/case03/input/input.txt index fbe8bcb..40b0a9d 100644 --- a/examples/case03/input/input.txt +++ b/examples/case03/input/input.txt @@ -8,6 +8,7 @@ step_simulate: 1 # 运行物理模拟 → output/trajectory.txt step_sample: 1 # 抽帧 → output/display.txt step_plot: 0 # 绘制轨迹/能量图 → output/trajectory_plots.png +step_plot_wave: 0 # 绘制波形能量动画 → output/wave_animation.gif step_animation: 1 # 自动播放 VisPy 3D 动画窗口(需安装 vispy) force_calc: 0 # 强制重新计算:1=跳过缓存强算,0=自动使用已有输出 diff --git a/examples/case04/input/input.txt b/examples/case04/input/input.txt index fbe8bcb..40b0a9d 100644 --- a/examples/case04/input/input.txt +++ b/examples/case04/input/input.txt @@ -8,6 +8,7 @@ step_simulate: 1 # 运行物理模拟 → output/trajectory.txt step_sample: 1 # 抽帧 → output/display.txt step_plot: 0 # 绘制轨迹/能量图 → output/trajectory_plots.png +step_plot_wave: 0 # 绘制波形能量动画 → output/wave_animation.gif step_animation: 1 # 自动播放 VisPy 3D 动画窗口(需安装 vispy) force_calc: 0 # 强制重新计算:1=跳过缓存强算,0=自动使用已有输出 diff --git a/examples/case05/doc/index.html b/examples/case05/doc/index.html index 07267ac..88e8b7d 100644 --- a/examples/case05/doc/index.html +++ b/examples/case05/doc/index.html @@ -356,10 +356,14 @@ step_animation: 1 # 播放动画 操作效果 鼠标拖动旋转视角 滚轮缩放 + W / S 键相机沿 Z 轴向前 / 向后移动(靠近/远离场景) + A / D 键视角向右 / 向左平移 + E / Q 键视角上升 / 下降(屏幕方向) + C / X 键增大 / 减小步长 + V 键切换透视 / 正交投影 左上角 reset 按钮复位视角到初始位置 左上角 info 按钮切换信息面板显示/隐藏 左上角 axes 按钮切换坐标轴显示/隐藏 - Q / E 键画面绕视向旋转 -90° / +90° @@ -395,6 +399,7 @@ step_animation: 1 # 播放动画 step_simulate跳过模拟(加载已有轨迹)运行物理模拟 step_sample跳过抽帧从轨迹抽取显示帧 step_plot不生成图表生成轨迹/能量图 + step_plot_wave不生成波形图生成波形能量动画 GIF step_animation不启动动画自动打开 VisPy 3D 窗口 force_calc自动检测缓存强制重新计算 @@ -418,7 +423,8 @@ step_animation: 1 # 播放动画 │ ├── trajectory.txt # 全量轨迹数据(10000 步 × 60 原子) │ ├── display.txt # 抽帧后的动画数据(200 帧 × 60 原子) │ ├── dynamics.log # 计算日志 -│ └── animation.log # 动画启动日志(闪退时排查用) +│ ├── animation.log # 动画启动日志(闪退时排查用) +│ └── wave_animation.gif # 波形能量动画(step_plot_wave=1 时生成) ├── doc/ │ └── index.html # 本文档 ├── Readme.md # 案例简介 diff --git a/examples/case05/input/bond.txt b/examples/case05/input/bond.txt index f33e17f..4434bf4 100644 --- a/examples/case05/input/bond.txt +++ b/examples/case05/input/bond.txt @@ -1,2 +1,2 @@ bond_name k rest_length -k1 50.0 1.0 +k1 100.0 1.0 diff --git a/examples/case05/input/coord.txt b/examples/case05/input/coord.txt index f34b022..8bd9dda 100644 --- a/examples/case05/input/coord.txt +++ b/examples/case05/input/coord.txt @@ -1,61 +1,61 @@ n mass radius x y z vx vy vz fix_x fix_y fix_z -1 1 0.1 0 0 1 0 0 0 1 1 0 -2 1 0.1 1 0 0 0 0 0 1 1 0 -3 1 0.1 2 0 0 0 0 0 1 1 0 -4 1 0.1 3 0 0 0 0 0 1 1 0 -5 1 0.1 4 0 0 0 0 0 1 1 0 -6 1 0.1 5 0 0 0 0 0 1 1 0 -7 1 0.1 6 0 0 0 0 0 1 1 0 -8 1 0.1 7 0 0 0 0 0 1 1 0 -9 1 0.1 8 0 0 0 0 0 1 1 0 -10 1 0.1 9 0 0 0 0 0 1 1 0 -11 1 0.1 10 0 0 0 0 0 1 1 0 -12 1 0.1 11 0 0 0 0 0 1 1 0 -13 1 0.1 12 0 0 0 0 0 1 1 0 -14 1 0.1 13 0 0 0 0 0 1 1 0 -15 1 0.1 14 0 0 0 0 0 1 1 0 -16 1 0.1 15 0 0 0 0 0 1 1 0 -17 1 0.1 16 0 0 0 0 0 1 1 0 -18 1 0.1 17 0 0 0 0 0 1 1 0 -19 1 0.1 18 0 0 0 0 0 1 1 0 -20 1 0.1 19 0 0 0 0 0 1 1 0 -21 1 0.1 20 0 0 0 0 0 1 1 0 -22 1 0.1 21 0 0 0 0 0 1 1 0 -23 1 0.1 22 0 0 0 0 0 1 1 0 -24 1 0.1 23 0 0 0 0 0 1 1 0 -25 1 0.1 24 0 0 0 0 0 1 1 0 -26 1 0.1 25 0 0 0 0 0 1 1 0 -27 1 0.1 26 0 0 0 0 0 1 1 0 -28 1 0.1 27 0 0 0 0 0 1 1 0 -29 1 0.1 28 0 0 0 0 0 1 1 0 -30 1 0.1 29 0 0 0 0 0 1 1 0 -31 1 0.1 30 0 0 0 0 0 1 1 0 -32 1 0.1 31 0 0 0 0 0 1 1 0 -33 1 0.1 32 0 0 0 0 0 1 1 0 -34 1 0.1 33 0 0 0 0 0 1 1 0 -35 1 0.1 34 0 0 0 0 0 1 1 0 -36 1 0.1 35 0 0 0 0 0 1 1 0 -37 1 0.1 36 0 0 0 0 0 1 1 0 -38 1 0.1 37 0 0 0 0 0 1 1 0 -39 1 0.1 38 0 0 0 0 0 1 1 0 -40 1 0.1 39 0 0 0 0 0 1 1 0 -41 1 0.1 40 0 0 0 0 0 1 1 0 -42 1 0.1 41 0 0 0 0 0 1 1 0 -43 1 0.1 42 0 0 0 0 0 1 1 0 -44 1 0.1 43 0 0 0 0 0 1 1 0 -45 1 0.1 44 0 0 0 0 0 1 1 0 -46 1 0.1 45 0 0 0 0 0 1 1 0 -47 1 0.1 46 0 0 0 0 0 1 1 0 -48 1 0.1 47 0 0 0 0 0 1 1 0 -49 1 0.1 48 0 0 0 0 0 1 1 0 -50 1 0.1 49 0 0 0 0 0 1 1 0 -51 1 0.1 50 0 0 0 0 0 1 1 0 -52 1 0.1 51 0 0 0 0 0 1 1 0 -53 1 0.1 52 0 0 0 0 0 1 1 0 -54 1 0.1 53 0 0 0 0 0 1 1 0 -55 1 0.1 54 0 0 0 0 0 1 1 0 -56 1 0.1 55 0 0 0 0 0 1 1 0 -57 1 0.1 56 0 0 0 0 0 1 1 0 -58 1 0.1 57 0 0 0 0 0 1 1 0 -59 1 0.1 58 0 0 0 0 0 1 1 0 -60 1 0.1 59 0 0 0 0 0 1 1 0 +1 1 0.1 0 0 1 0 0 0 0 1 1 +2 1 0.1 1 0 0 0 0 0 0 1 1 +3 1 0.1 2 0 0 0 0 0 0 1 1 +4 1 0.1 3 0 0 0 0 0 0 1 1 +5 1 0.1 4 0 0 0 0 0 0 1 1 +6 1 0.1 5 0 0 0 0 0 0 1 1 +7 1 0.1 6 0 0 0 0 0 0 1 1 +8 1 0.1 7 0 0 0 0 0 0 1 1 +9 1 0.1 8 0 0 0 0 0 0 1 1 +10 1 0.1 9 0 0 0 0 0 0 1 1 +11 1 0.1 10 0 0 0 0 0 0 1 1 +12 1 0.1 11 0 0 0 0 0 0 1 1 +13 1 0.1 12 0 0 0 0 0 0 1 1 +14 1 0.1 13 0 0 0 0 0 0 1 1 +15 1 0.1 14 0 0 0 0 0 0 1 1 +16 1 0.1 15 0 0 0 0 0 0 1 1 +17 1 0.1 16 0 0 0 0 0 0 1 1 +18 1 0.1 17 0 0 0 0 0 0 1 1 +19 1 0.1 18 0 0 0 0 0 0 1 1 +20 1 0.1 19 0 0 0 0 0 0 1 1 +21 1 0.1 20 0 0 0 0 0 0 1 1 +22 1 0.1 21 0 0 0 0 0 0 1 1 +23 1 0.1 22 0 0 0 0 0 0 1 1 +24 1 0.1 23 0 0 0 0 0 0 1 1 +25 1 0.1 24 0 0 0 0 0 0 1 1 +26 1 0.1 25 0 0 0 0 0 0 1 1 +27 1 0.1 26 0 0 0 0 0 0 1 1 +28 1 0.1 27 0 0 0 0 0 0 1 1 +29 1 0.1 28 0 0 0 0 0 0 1 1 +30 1 0.1 29 0 0 0 0 0 0 1 1 +31 1 0.1 30 0 0 0 0 0 0 1 1 +32 1 0.1 31 0 0 0 0 0 0 1 1 +33 1 0.1 32 0 0 0 0 0 0 1 1 +34 1 0.1 33 0 0 0 0 0 0 1 1 +35 1 0.1 34 0 0 0 0 0 0 1 1 +36 1 0.1 35 0 0 0 0 0 0 1 1 +37 1 0.1 36 0 0 0 0 0 0 1 1 +38 1 0.1 37 0 0 0 0 0 0 1 1 +39 1 0.1 38 0 0 0 0 0 0 1 1 +40 1 0.1 39 0 0 0 0 0 0 1 1 +41 1 0.1 40 0 0 0 0 0 0 1 1 +42 1 0.1 41 0 0 0 0 0 0 1 1 +43 1 0.1 42 0 0 0 0 0 0 1 1 +44 1 0.1 43 0 0 0 0 0 0 1 1 +45 1 0.1 44 0 0 0 0 0 0 1 1 +46 1 0.1 45 0 0 0 0 0 0 1 1 +47 1 0.1 46 0 0 0 0 0 0 1 1 +48 1 0.1 47 0 0 0 0 0 0 1 1 +49 1 0.1 48 0 0 0 0 0 0 1 1 +50 1 0.1 49 0 0 0 0 0 0 1 1 +51 1 0.1 50 0 0 0 0 0 0 1 1 +52 1 0.1 51 0 0 0 0 0 0 1 1 +53 1 0.1 52 0 0 0 0 0 0 1 1 +54 1 0.1 53 0 0 0 0 0 0 1 1 +55 1 0.1 54 0 0 0 0 0 0 1 1 +56 1 0.1 55 0 0 0 0 0 0 1 1 +57 1 0.1 56 0 0 0 0 0 0 1 1 +58 1 0.1 57 0 0 0 0 0 0 1 1 +59 1 0.1 58 0 0 0 0 0 0 1 1 +60 1 0.1 59 0 0 0 0 0 1 1 1 diff --git a/examples/case05/input/driver.txt b/examples/case05/input/driver.txt index 79857fc..2744c96 100644 --- a/examples/case05/input/driver.txt +++ b/examples/case05/input/driver.txt @@ -1,2 +1,2 @@ n amp_x amp_y amp_z freq_x freq_y freq_z phi_x phi_y phi_z period -1 0 0 1 0 0 0.1 0 0 90 all \ No newline at end of file +1 0.5 0 0 1 0 0 0 0 0 all \ No newline at end of file diff --git a/examples/case05/input/input.txt b/examples/case05/input/input.txt index 039356c..38afc3b 100644 --- a/examples/case05/input/input.txt +++ b/examples/case05/input/input.txt @@ -5,15 +5,16 @@ # ── 流程控制 ────────────────────────────────── # 每步用 0/1 单独开关,1=执行,0=跳过 # 依赖关系:抽帧依赖模拟结果,绘图依赖模拟+抽帧 -step_simulate: 1 # 运行物理模拟 → output/trajectory.txt -step_sample: 1 # 抽帧 → output/display.txt +step_simulate: 0 # 运行物理模拟 → output/trajectory.txt +step_sample: 0 # 抽帧 → output/display.txt step_plot: 0 # 绘制轨迹/能量图 → output/trajectory_plots.png -step_animation: 1 # 自动播放 VisPy 3D 动画窗口(需安装 vispy) +step_animation: 0 # 自动播放 VisPy 3D 动画窗口(需安装 vispy) +step_plot_wave: 1 # 绘制波形图 force_calc: 0 # 强制重新计算:1=跳过缓存强算,0=自动使用已有输出 # ── 计算引擎 ────────────────────────────────── # 可选: python, c, cpp, fortran, java -engine: python # 默认使用 Python 引擎 +engine: python # 默认使用 python 引擎 # ── 盒子 ────────────────────────────────────── box_a: 80.0 # 立方体半边长,粒子被限制在 [-box_a, box_a]³ 内 @@ -60,7 +61,7 @@ warmup_steps: 0 # 默认 0(立即开始记录) # 总模拟时间(秒),程序自动计算 NT = T_total / DT # 如果同时指定了 NT,以 NT 为准 -T_total: 100.0 +T_total: 10.0 # 抽帧间隔(每 NSTEP 步取一帧用于动画) NSTEP: 50 diff --git a/examples/case06/Readme.md b/examples/case06/Readme.md new file mode 100644 index 0000000..2422398 --- /dev/null +++ b/examples/case06/Readme.md @@ -0,0 +1,42 @@ +# case05: 一维原子链驱动力学模拟 + +60 个原子沿 x 轴排列,相邻原子用弹簧连接。原子 1 受驱动力作用。 + +## 物理设定 + +| 参数 | 值 | +|---|---| +| 原子数 | 60 | +| 排列 | 沿 x 轴等间距排列,间距为 1 | +| 约束 | 原子**只沿 z 方向**振动(fix_x=1, fix_y=1, fix_z=0) | +| 弹簧 | 劲度系数 k=1.0,原长 L₀=1.0 | +| 重力 | 无 | +| 万有引力 | 无 | +| 阻尼 | 无 | +| 驱动力 | 原子 1(详见 driver.txt) | +| 算法 | leapfrog(蛙跳法,能量守恒) | + +## 驱动力 + +原子 1 的位置由 `input/driver.txt` 中的驱动力公式决定: + +```math +z(t) = A_z \cdot \cos(2\pi f_z t + \phi_z) +``` + +当前参数:A_z = 5.0, f_z = 1.0 Hz, φ_z = 90°(全程驱动)。 + +受驱原子完全忽略 coord.txt 中的初始坐标和 fix 约束,位置/速度由驱动力解析确定。 + +## 动力学行为 + +原子 1 的受迫振动通过弹簧逐次传递给相邻原子,形成沿链传播的横波。由于横向振动的几何非线性(弹簧张力主要在 x 方向),z 方向有效刚度低,波速较慢,呈现 FPU 型非线性动力学特征。 + +## 使用方法 + +```bash +cd examples/case05 +python run_dynamics.py +``` + +配置参数详见 `input/input.txt`,驱动力定义见 `input/driver.txt`,完整文档见 `doc/index.html`。 diff --git a/examples/case06/doc/index.html b/examples/case06/doc/index.html new file mode 100644 index 0000000..7aa9e48 --- /dev/null +++ b/examples/case06/doc/index.html @@ -0,0 +1,473 @@ + + + + + +case05 — 一维原子链驱动力学模拟 | 物理原理 & 使用文档 + + + + + + + +
+

一维原子链驱动力学模拟

+

60 个原子沿 x 轴排列 · 弹簧连接 · z 方向受迫振动

+ case05 · examples/case05 +
+ +
+ + + + +
+

目录

+
    +
  1. 物理原理
  2. +
  3. 数值算法
  4. +
  5. 驱动力模型
  6. +
  7. 使用方法
  8. +
  9. 参数参考
  10. +
  11. 文件结构
  12. +
  13. 常见问题
  14. +
+
+ + + + +
+

一、物理原理

+ +
+

1.1 一维原子链

+

60 个原子沿 x 轴 等间距排列,原子间距为 1。相邻原子之间用 理想弹簧 连接,弹簧的劲度系数 k = 1.0,原长 L₀ = 1.0(与原子间距一致,初始状态弹簧无拉伸)。

+

每个原子被限制在 z 方向 自由振动,x 和 y 方向锁定(fix_x=1, fix_y=1, fix_z=0)。

+
+ +
+

1.2 弹簧力(胡克定律)

+

当原子 ij 之间有弹簧连接时,原子 i 受到的弹簧力为:

+
+ F = −k · (dL₀) · uij +
+

其中 d = |rjri| 为两原子间距离,uij 为从 i 指向 j 的单位向量。由于原子只在 z 方向振动,弹簧在 z 方向的分量是 几何非线性 的——对于小振幅近似,z 方向等效于一个三次方恢复力(FPU 型非线性)。

+
+ +
+

1.3 运动方程

+

对于第 i 个自由原子(非受驱),牛顿第二定律给出:

+
+ m · ai = Fispring + Fidriving +
+

本案例中 唯一的外力 来自驱动力(仅施加于原子 1)。无重力、无万有引力、无阻尼,系统总能量守恒。

+
+ +
+

1.4 波传播

+

原子 1 的受迫振动通过弹簧逐次传递给相邻原子,形成沿链传播的 横波。由于横向振动的几何非线性(弹簧大部分张力在 x 方向,z 方向的有效刚度远小于 1),波的传播速度较慢,且高阶频率成分会在链中产生复杂的非线性动力学行为(类似 FPU 回波现象)。

+
+
+ + + + +
+

二、数值算法

+ +
+

2.1 蛙跳法(Leapfrog / Velocity-Verlet)

+

采用能量守恒特性优异的 蛙跳法(二阶辛积分器),更新公式为:

+
+ v(t + ½Δt) = v(t) + ½ a(t) · Δt
+ r(t + Δt) = r(t) + v(t + ½Δt) · Δt
+ a(t + Δt) = F(r(t + Δt), v(t + ½Δt)) / m
+ v(t + Δt) = v(t + ½Δt) + ½ a(t + Δt) · Δt +
+

蛙跳法在长时间模拟中能量漂移极小(本案例验证 < 0.004%),适合无阻尼的保守系统。

+
+ +
+

2.2 时间步长与采样

+ + + + + + +
参数说明
DT0.01 s积分步长(远小于 1/ω ≈ 0.16 s,满足稳定性条件)
T_total100 s总模拟时间 → NT = 10000 步
NSTEP50每 NSTEP 步取一帧用于动画 → 200 帧
methodleapfrog蛙跳法(Velocity-Verlet)
+
+ +
+

2.3 计算流程

+
+ 读入 coord.txt
connection.txt
bond.txt
+ + 施加驱动力
(驱动原子 1)
+ + 记录轨迹 + + 蛙跳法
更新位置/速度
+ + 固定约束
(x, y 锁定)
+ + 循环
NT 次
+
+

注意:驱动力在 每次积分前 施加,确保受驱原子的位置正确传递给弹簧力计算。

+
+
+ + + + +
+

三、驱动力模型

+ +
+

3.1 定义文件

+

驱动力由 input/driver.txt 定义,格式如下:

+
n amp_x amp_y amp_z freq_x freq_y freq_z phi_x phi_y phi_z period
+1     0     0     5      0      0      1     0     0    90    all
+
+ +
+

3.2 数学公式

+

受驱原子的位置由下式决定(完全替换 coord.txt 中的初始坐标和固定约束):

+
+ r(t) = A · cos(2πf · t + φ) +
+

速度由解析导数给出:

+
+ v(t) = −A · 2πf · sin(2πf · t + φ) +
+

其中 A = (amp_x, amp_y, amp_z),f = (freq_x, freq_y, freq_z) 为不同方向的驱动频率,φ = (phi_x, phi_y, phi_z) 为相位(角度制,代码自动转换为弧度)。

+
+ +
+

3.3 本案例驱动参数

+ + + + + + +
参数含义
amp_z5.0z 方向驱动振幅
freq_z1.0 Hz驱动频率(周期 1 s)
phi_z90°驱动相位 → z(0) = 5·cos(90°) = 0
periodall全程驱动,永不停止
+
+ z(t) = 5.0 · cos(2π · 1.0 · t + 90°) +
+
+ +
+

3.4 有限周期驱动

+

period 参数支持三种模式:

+
    +
  • all — 全程驱动
  • +
  • 数值 — 驱动指定周期数后 静止(冻结在最终位置,速度归零)。例如 period: 1 表示驱动 1 个完整周期后停止。
  • +
+
+ +
+

3.5 驱动与固定约束的关系

+

对于受驱原子(driver.txtn 指定的原子),其在 coord.txt 中的初始坐标和 fix_x/fix_y/fix_z 约束被 完全忽略。原子的位置和速度完全由驱动力公式决定。

+
+
+ + + + +
+

四、使用方法

+ +
+

4.1 完整运行(模拟 + 动画)

+
cd examples/case05
+python run_dynamics.py
+

这步会依次执行:物理模拟 → 抽帧 → 打开 VisPy 3D 动画窗口。

+
+ +
+

4.2 仅查看已有结果

+

如果已经跑完模拟且生成了 output/display.txt,可以通过修改 input.txt 跳过计算,只开动画:

+
step_simulate:  0   # 跳过模拟
+step_sample:    0   # 跳过抽帧
+step_animation: 1   # 播放动画
+

然后运行:python run_dynamics.py

+
+ +
+

4.3 手动 3D 动画

+

也可以单独启动 VisPy 窗口:

+
python ../../draw.py output/
+
+ +
+

4.4 强制重新计算

+

修改参数后需要重新运行模拟时,设置:

+
force_calc: 1        # 忽略缓存,强制重新计算
+
+ +
+

4.5 动画交互

+ + + + + + + + + + +
操作效果
鼠标拖动旋转视角
滚轮缩放
A / D 键视角向左 / 向右平移
W / S 键视角向上 / 向下平移
Q / E 键降低 / 提高平移速度
左上角 reset 按钮复位视角到初始位置
左上角 info 按钮切换信息面板显示/隐藏
左上角 axes 按钮切换坐标轴显示/隐藏
+
+
+ + + + +
+

五、参数参考

+ +
+

5.1 input.txt 关键参数

+ + + + + + + + + + + + + +
参数默认值说明
gravity_field0均匀重力场(已关闭)
gravity_interaction0原子间万有引力(已关闭)
elastic_force1弹簧键力(已开启)
damping_force0阻尼(已关闭)
driving_force1驱动力开关(1=开启,需 driver.txt)
methodleapfrog数值积分方法
DT0.01积分步长 (s)
T_total100.0总模拟时间 (s)
NSTEP50抽帧步数间隔
enginepython计算引擎(python / c / cpp / fortran)
use_marker1渲染模式(0=Sphere 网格, 1=Marker GPU 实例化)
+
+ +
+

5.2 流程控制参数

+ + + + + + + +
参数01
step_simulate跳过模拟(加载已有轨迹)运行物理模拟
step_sample跳过抽帧从轨迹抽取显示帧
step_plot不生成图表生成轨迹/能量图
step_animation不启动动画自动打开 VisPy 3D 窗口
force_calc自动检测缓存强制重新计算
+
+
+ + + + +
+

六、文件结构

+ +
case05/
+├── input/
+│   ├── input.txt          # 主配置文件(YAML 格式)
+│   ├── coord.txt          # 原子坐标(60 个原子)
+│   ├── connection.txt     # 弹簧连接关系(59 条键)
+│   ├── bond.txt           # 弹簧参数(k=1.0, L₀=1.0)
+│   └── driver.txt       # 驱动力定义(本案例新增)
+├── output/
+│   ├── trajectory.txt     # 全量轨迹数据(10000 步 × 60 原子)
+│   ├── display.txt        # 抽帧后的动画数据(200 帧 × 60 原子)
+│   ├── dynamics.log       # 计算日志
+│   └── animation.log      # 动画启动日志(闪退时排查用)
+├── doc/
+│   └── index.html         # 本文档
+├── Readme.md              # 案例简介
+└── run_dynamics.py        # 案例运行入口
+
+ + + + +
+

七、常见问题

+ +
+

7.1 动画窗口闪退

+

如果 VisPy 窗口一闪就消失,请检查:

+
    +
  • output/animation.log 中是否有错误信息
  • +
  • output/display.txt 是否存在(需先跑 step_sample: 1
  • +
+
+ +
+

7.2 原子不振动

+

可能原因:

+
    +
  • NSTEP 过大:抽帧间隔大于驱动周期的一半时,动画会丢失振动细节。建议 NSTEP ≤ 1/(freq · DT · 10)
  • +
  • 相位 φ 使采样点落在零值:试试 phi_z: 0 让原子在 t=0 处于振幅峰值
  • +
  • 确认 driving_force: 1driver.txt 中 amp_z 不为 0
  • +
+
+ +
+

7.3 渲染性能慢

+

原子数多时动画卡顿:

+
    +
  • 设置 use_marker: 1(使用 GPU 实例化渲染替代独立网格球体)
  • +
  • 增大 NSTEP 减少动画帧数
  • +
+
+
+ +
+ + + +
+ + diff --git a/examples/case06/input/bond.txt b/examples/case06/input/bond.txt new file mode 100644 index 0000000..4434bf4 --- /dev/null +++ b/examples/case06/input/bond.txt @@ -0,0 +1,2 @@ +bond_name k rest_length +k1 100.0 1.0 diff --git a/examples/case06/input/connection.txt b/examples/case06/input/connection.txt new file mode 100644 index 0000000..c0d1094 --- /dev/null +++ b/examples/case06/input/connection.txt @@ -0,0 +1,60 @@ +n1 n2 bond_name +1 2 k1 +2 3 k1 +3 4 k1 +4 5 k1 +5 6 k1 +6 7 k1 +7 8 k1 +8 9 k1 +9 10 k1 +10 11 k1 +11 12 k1 +12 13 k1 +13 14 k1 +14 15 k1 +15 16 k1 +16 17 k1 +17 18 k1 +18 19 k1 +19 20 k1 +20 21 k1 +21 22 k1 +22 23 k1 +23 24 k1 +24 25 k1 +25 26 k1 +26 27 k1 +27 28 k1 +28 29 k1 +29 30 k1 +30 31 k1 +31 32 k1 +32 33 k1 +33 34 k1 +34 35 k1 +35 36 k1 +36 37 k1 +37 38 k1 +38 39 k1 +39 40 k1 +40 41 k1 +41 42 k1 +42 43 k1 +43 44 k1 +44 45 k1 +45 46 k1 +46 47 k1 +47 48 k1 +48 49 k1 +49 50 k1 +50 51 k1 +51 52 k1 +52 53 k1 +53 54 k1 +54 55 k1 +55 56 k1 +56 57 k1 +57 58 k1 +58 59 k1 +59 60 k1 diff --git a/examples/case06/input/coord.txt b/examples/case06/input/coord.txt new file mode 100644 index 0000000..8bd9dda --- /dev/null +++ b/examples/case06/input/coord.txt @@ -0,0 +1,61 @@ +n mass radius x y z vx vy vz fix_x fix_y fix_z +1 1 0.1 0 0 1 0 0 0 0 1 1 +2 1 0.1 1 0 0 0 0 0 0 1 1 +3 1 0.1 2 0 0 0 0 0 0 1 1 +4 1 0.1 3 0 0 0 0 0 0 1 1 +5 1 0.1 4 0 0 0 0 0 0 1 1 +6 1 0.1 5 0 0 0 0 0 0 1 1 +7 1 0.1 6 0 0 0 0 0 0 1 1 +8 1 0.1 7 0 0 0 0 0 0 1 1 +9 1 0.1 8 0 0 0 0 0 0 1 1 +10 1 0.1 9 0 0 0 0 0 0 1 1 +11 1 0.1 10 0 0 0 0 0 0 1 1 +12 1 0.1 11 0 0 0 0 0 0 1 1 +13 1 0.1 12 0 0 0 0 0 0 1 1 +14 1 0.1 13 0 0 0 0 0 0 1 1 +15 1 0.1 14 0 0 0 0 0 0 1 1 +16 1 0.1 15 0 0 0 0 0 0 1 1 +17 1 0.1 16 0 0 0 0 0 0 1 1 +18 1 0.1 17 0 0 0 0 0 0 1 1 +19 1 0.1 18 0 0 0 0 0 0 1 1 +20 1 0.1 19 0 0 0 0 0 0 1 1 +21 1 0.1 20 0 0 0 0 0 0 1 1 +22 1 0.1 21 0 0 0 0 0 0 1 1 +23 1 0.1 22 0 0 0 0 0 0 1 1 +24 1 0.1 23 0 0 0 0 0 0 1 1 +25 1 0.1 24 0 0 0 0 0 0 1 1 +26 1 0.1 25 0 0 0 0 0 0 1 1 +27 1 0.1 26 0 0 0 0 0 0 1 1 +28 1 0.1 27 0 0 0 0 0 0 1 1 +29 1 0.1 28 0 0 0 0 0 0 1 1 +30 1 0.1 29 0 0 0 0 0 0 1 1 +31 1 0.1 30 0 0 0 0 0 0 1 1 +32 1 0.1 31 0 0 0 0 0 0 1 1 +33 1 0.1 32 0 0 0 0 0 0 1 1 +34 1 0.1 33 0 0 0 0 0 0 1 1 +35 1 0.1 34 0 0 0 0 0 0 1 1 +36 1 0.1 35 0 0 0 0 0 0 1 1 +37 1 0.1 36 0 0 0 0 0 0 1 1 +38 1 0.1 37 0 0 0 0 0 0 1 1 +39 1 0.1 38 0 0 0 0 0 0 1 1 +40 1 0.1 39 0 0 0 0 0 0 1 1 +41 1 0.1 40 0 0 0 0 0 0 1 1 +42 1 0.1 41 0 0 0 0 0 0 1 1 +43 1 0.1 42 0 0 0 0 0 0 1 1 +44 1 0.1 43 0 0 0 0 0 0 1 1 +45 1 0.1 44 0 0 0 0 0 0 1 1 +46 1 0.1 45 0 0 0 0 0 0 1 1 +47 1 0.1 46 0 0 0 0 0 0 1 1 +48 1 0.1 47 0 0 0 0 0 0 1 1 +49 1 0.1 48 0 0 0 0 0 0 1 1 +50 1 0.1 49 0 0 0 0 0 0 1 1 +51 1 0.1 50 0 0 0 0 0 0 1 1 +52 1 0.1 51 0 0 0 0 0 0 1 1 +53 1 0.1 52 0 0 0 0 0 0 1 1 +54 1 0.1 53 0 0 0 0 0 0 1 1 +55 1 0.1 54 0 0 0 0 0 0 1 1 +56 1 0.1 55 0 0 0 0 0 0 1 1 +57 1 0.1 56 0 0 0 0 0 0 1 1 +58 1 0.1 57 0 0 0 0 0 0 1 1 +59 1 0.1 58 0 0 0 0 0 0 1 1 +60 1 0.1 59 0 0 0 0 0 1 1 1 diff --git a/examples/case06/input/driver.txt b/examples/case06/input/driver.txt new file mode 100644 index 0000000..69269f3 --- /dev/null +++ b/examples/case06/input/driver.txt @@ -0,0 +1,2 @@ +n amp_x amp_y amp_z freq_x freq_y freq_z phi_x phi_y phi_z period +1 0.5 0 0 0.1 0 0 0 0 0 all \ No newline at end of file diff --git a/examples/case06/input/input.txt b/examples/case06/input/input.txt new file mode 100644 index 0000000..d52a9fa --- /dev/null +++ b/examples/case06/input/input.txt @@ -0,0 +1,96 @@ +# 物理模拟参数配置 +# 格式:YAML +# 用法:python run_dynamics.py + +# ── 流程控制 ────────────────────────────────── +# 每步用 0/1 单独开关,1=执行,0=跳过 +# 依赖关系:抽帧依赖模拟结果,绘图依赖模拟+抽帧 +step_simulate: 1 # 运行物理模拟 → output/trajectory.txt +step_sample: 1 # 抽帧 → output/display.txt +step_plot: 0 # 绘制轨迹/能量图 → output/trajectory_plots.png +step_animation: 1 # 自动播放 VisPy 3D 动画窗口(需安装 vispy) +force_calc: 0 # 强制重新计算:1=跳过缓存强算,0=自动使用已有输出 + +# ── 计算引擎 ────────────────────────────────── +# 可选: python, c, cpp, fortran, java +engine: python # 默认使用 python 引擎 + +# ── 盒子 ────────────────────────────────────── +box_a: 80.0 # 立方体半边长,粒子被限制在 [-box_a, box_a]³ 内 + +# ── 初始构型 ────────────────────────────────── +# 坐标文件格式: +# 第一行:n mass radius x y z vx vy vz fix_x fix_y fix_z +# 后续行:原子序号 质量 半径 x y z vx vy vz fix_x fix_y fix_z +coord_file: input/coord.txt +connection_file: input/connection.txt +bond_file: input/bond.txt +driver_file: input/driver.txt # 驱动力定义文件(driving_force=1 时生效) + +# 绘图/动画展示的原子序号(对应 coord_file 第一列 n) +plot_atom: 1 + +# ── 物理参数 ────────────────────────────────── +# 三个方向分量分别对应 x, y, z +G: [0.0, 0.0, 0.0] # 重力场分量 (m/s²) +B: [0.0, 0.0, 0.0] # 阻尼分量 + +# ── 力开关(0=关闭, 1=开启)────────────────── +gravity_field: 0 # 均匀重力场 (G) +gravity_interaction: 0 # 原子间万有引力 +elastic_force: 1 # 弹簧键力 +damping_force: 0 # 阻尼 (B) +driving_force: 1 # 驱动力(需 driver_file 定义) +# +gravity_strength: 1.0 # 万有引力强度(仅 gravity_interaction=1 时有效) + +# ── 数值算法 ────────────────────────────────── +# 可选: +# explicit_euler 显式欧拉法 +# implicit_euler 隐式欧拉法 +# midpoint 中点法 +# leapfrog 蛙跳法 +method: leapfrog + +# ── 步骤控制 ────────────────────────────────── +# 以下参数控制哪些步骤被执行和保存 + +# 预热步数:模拟开始时跳过不保存的步数(用于稳定初始状态) +warmup_steps: 0 # 默认 0(立即开始记录) + +# 总模拟时间(秒),程序自动计算 NT = T_total / DT +# 如果同时指定了 NT,以 NT 为准 +T_total: 50.0 + +# 抽帧间隔(每 NSTEP 步取一帧用于动画) +NSTEP: 50 + +# ── 时间步长 ────────────────────────────────── +DT: 0.001 # 时间步长 (s) + +# 抽帧范围:只保存 [sample_start, sample_end) 区间内的帧 +sample_start: null # null 表示从头开始(帧索引从 0 起) +sample_end: null # null 表示到末尾 + + + +# ── 渲染方式 ────────────────────────────────── +# 3D 动画中原子渲染方式: +# 0 = Sphere (网格球体,效果精细,原子数少时推荐) +# 1 = Marker (GPU 实例化点,原子数多时性能更佳) +use_marker: 1 + +# ── 显示参数 ────────────────────────────────── +# 盒子透明度:单个数值(统一)或 6 个数的数组,按 [-x,+x,-y,+y,-z,+z] 顺序 +alpha: [0.0, 0.0, 0.0, 0.0, 0.0, 0.5] + +# 小球颜色 +# 小球半径从 coord_file 的 radius 列读取 +ball_color_r: 0.20 # R 分量 (0~1) +ball_color_g: 0.60 # G 分量 +ball_color_b: 0.90 # B 分量 + +# 盒子面颜色 +box_color_r: 0.80 +box_color_g: 0.80 +box_color_b: 0.85 diff --git a/examples/case06/run_dynamics.py b/examples/case06/run_dynamics.py new file mode 100644 index 0000000..37c77d2 --- /dev/null +++ b/examples/case06/run_dynamics.py @@ -0,0 +1,54 @@ +""" +Case runner for Dynamics case05 — 1D atomic chain. + +This script keeps program and data separated: + - program: ../../dynamics.py + - input: ./input + - output: ./output +""" + +from __future__ import annotations + +import argparse +import importlib.util +from pathlib import Path + + +CASE_DIR = Path(__file__).resolve().parent +DYNAMICS_PATH = Path("..") / ".." / "dynamics.py" +INPUT_DIR = Path("input") +OUTPUT_DIR = Path("output") +CONFIG_FILE = INPUT_DIR / "input.txt" + + +def load_dynamics_module(module_path: Path): + spec = importlib.util.spec_from_file_location("dynamics_module", module_path) + if spec is None or spec.loader is None: + raise ImportError(f"无法加载 dynamics.py: {module_path}") + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + return module + + +def main(): + parser = argparse.ArgumentParser(description="运行 Dynamics 示例案例 case01") + parser.add_argument("--no-plot", action="store_true", help="跳过 matplotlib 绘图") + args = parser.parse_args() + + dynamics_path = (CASE_DIR / DYNAMICS_PATH).resolve() + input_dir = (CASE_DIR / INPUT_DIR).resolve() + output_dir = (CASE_DIR / OUTPUT_DIR).resolve() + config_path = (CASE_DIR / CONFIG_FILE).resolve() + + module = load_dynamics_module(dynamics_path) + module.run_case( + config_path=config_path, + runtime_base=CASE_DIR, + input_dir=input_dir, + output_dir=output_dir, + no_plot=args.no_plot, + ) + + +if __name__ == "__main__": + main() diff --git a/plot_wave.py b/plot_wave.py new file mode 100644 index 0000000..0020318 --- /dev/null +++ b/plot_wave.py @@ -0,0 +1,260 @@ +""" +plot_wave.py +============ +波形与能量动态图:读取 display.txt,绘制原子位移波形 +(纵波 + 2 个横波)和系统能量/输入功率随时间变化的二维动画。 + +用法: + python plot_wave.py # 使用 dynamics 根目录下 output/ + python plot_wave.py examples/case05/output # 指定案例输出目录 +""" + +import numpy as np +import matplotlib.pyplot as plt +from matplotlib.animation import FuncAnimation +import os +import sys + +sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) +import compute + + +def load_disp_data(output_dir): + """加载 display.txt""" + disp_path = os.path.join(output_dir, "display.txt") + if not os.path.exists(disp_path): + raise FileNotFoundError(f"找不到 {disp_path}") + return compute.load_text_data(disp_path) + + +def compute_energy(x, y, z, vx, vy, vz, masses, mass_arr, + bond_pairs, bond_stiffness, bond_rest_lengths, + gravity_field, G, gravity_interaction, gravity_strength): + """计算系统各能量分量。 + + Returns: + ek_sys: 系统动能 (n_frames,) + us_sys: 系统弹性势能 (n_frames,) + ug_sys: 系统重力势能 (n_frames,) + ugr_sys: 系统万有引力势能 (n_frames,) + """ + n_frames = x.shape[0] + masses_2d = masses[np.newaxis, :] # (1, n_atoms) + + # 动能 Ek = ½ m v² + ek = 0.5 * masses_2d * (vx**2 + vy**2 + vz**2) + ek_sys = np.sum(ek, axis=1) + + # 弹性势能 Us = ½ k (d - d₀)² + us_sys = np.zeros(n_frames) + if bond_pairs is not None and len(bond_pairs) > 0: + for b_idx in range(len(bond_pairs)): + i, j = bond_pairs[b_idx] + dx = x[:, j] - x[:, i] + dy = y[:, j] - y[:, i] + dz = z[:, j] - z[:, i] + dist = np.sqrt(dx**2 + dy**2 + dz**2) + stretch = dist - bond_rest_lengths[b_idx] + us_sys += 0.5 * bond_stiffness[b_idx] * stretch**2 + + # 均匀重力场势能 Ug = -m G·r + ug_sys = np.zeros(n_frames) + if gravity_field: + G_vec = np.array(G) + ug_sys = -masses_2d * (G_vec[0] * x + G_vec[1] * y + G_vec[2] * z) + ug_sys = np.sum(ug_sys, axis=1) + + # 万有引力势能 Ug_grav = -G_grav Σ m_i m_j / r + ugr_sys = np.zeros(n_frames) + if gravity_interaction: + n_atoms = len(masses) + # 为避免巨大计算量,仅当原子数较少时计算 + if n_atoms <= 200: + for i in range(n_atoms): + for j in range(i + 1, n_atoms): + dx = x[:, j] - x[:, i] + dy = y[:, j] - y[:, i] + dz = z[:, j] - z[:, i] + dist = np.sqrt(dx**2 + dy**2 + dz**2) + dist = np.maximum(dist, 1e-12) + pair_pe = -gravity_strength * masses[i] * masses[j] / dist + ugr_sys += pair_pe + + return ek_sys, us_sys, ug_sys, ugr_sys + + +def plot_wave(output_dir): + """主绘图函数:读取 display.txt 并生成波形+能量动画。""" + data = load_disp_data(output_dir) + + n_frames = int(data["n_frames"]) + t = np.array(data["disp_t"]) + + # 位置 / 速度 + x = np.array(data["disp_all_x"]) + y = np.array(data["disp_all_y"]) + z = np.array(data["disp_all_z"]) + vx = np.array(data["disp_all_vx"]) + vy = np.array(data["disp_all_vy"]) + vz = np.array(data["disp_all_vz"]) + + # 原子信息 + pos_0 = np.array(data["atom_positions"]) # (n_atoms, 3) + masses = np.array(data["atom_masses"]) + atom_ids = np.array(data["atom_ids"]) + n_atoms = len(atom_ids) + + # 成键 + bond_pairs = data.get("bond_pairs", []) + bond_stiffness = np.array(data.get("bond_stiffness", [])) + bond_rest_lengths = np.array(data.get("bond_rest_lengths", [])) + + # 物理开关 + gravity_field = int(data.get("gravity_field", 0)) + gravity_interaction = int(data.get("gravity_interaction", 0)) + G = data.get("G", [0, 0, 0]) + gravity_strength = float(data.get("gravity_strength", 1.0)) + driving_force = int(data.get("driving_force", 0)) + + # ── 位移(偏离初始平衡位形)── + dx = x - pos_0[np.newaxis, :, 0] # 纵波(沿链方向 x) + dy = y - pos_0[np.newaxis, :, 1] # 横波 1(y 方向) + dz = z - pos_0[np.newaxis, :, 2] # 横波 2(z 方向) + + # ── 能量 ── + ek, us, ug, ugr = compute_energy( + x, y, z, vx, vy, vz, masses, masses, + bond_pairs, bond_stiffness, bond_rest_lengths, + gravity_field, G, gravity_interaction, gravity_strength) + e_total = ek + us + ug + ugr + power = np.gradient(e_total, t) # 输入功率 = dE/dt + + # ── 波形纵轴范围(全局统一,避免抖动)── + def get_ylim(disp_data): + vmax = np.max(np.abs(disp_data)) + if vmax < 1e-10: + return -1.0, 1.0 + margin = vmax * 0.2 + return -vmax - margin, vmax + margin + + ylims = [get_ylim(dx), get_ylim(dy), get_ylim(dz)] + e_max = max(np.max(e_total), 0.01) * 1.3 + p_abs = np.max(np.abs(power)) + p_max = max(p_abs * 1.3, 0.01) + + atom_idx = np.arange(n_atoms) + + # ── 图形设置 ── + plt.rcParams['font.sans-serif'] = ['Microsoft YaHei', 'SimHei', 'DejaVu Sans'] + plt.rcParams['axes.unicode_minus'] = False + + fig, axes = plt.subplots(2, 2, figsize=(14, 10)) + fig.suptitle("波形与能量分析", fontsize=16) + + # ── 3 个波形子图 ── + titles = [ + "纵波 (x 方向位移)", + "横波 (y 方向位移)", + "横波 (z 方向位移)", + ] + colors = ["#2563eb", "#ea580c", "#16a34a"] + wave_configs = list(zip( + [axes[0, 0], axes[0, 1], axes[1, 0]], + [dx, dy, dz], titles, colors, ylims + )) + + wave_lines = [] + time_texts = [] + for ax, disp, title, color, yl in wave_configs: + ax.set_xlim(0, n_atoms - 1) + ax.set_ylim(yl) + ax.set_xlabel("原子序号") + ax.set_ylabel("位移") + ax.set_title(title) + ax.grid(True, alpha=0.3) + line, = ax.plot([], [], color=color, linewidth=1.5) + wave_lines.append(line) + tt = ax.text(0.02, 0.95, "", transform=ax.transAxes, + fontsize=11, verticalalignment="top") + time_texts.append(tt) + + # ── 能量+功率子图 ── + ax_ep = axes[1, 1] + ax_ep.set_xlim(t[0], t[-1]) + ep_ylow = min(-0.1 * max(e_max, p_max), -p_max * 0.1) + ep_yhigh = max(e_max, p_max) + ax_ep.set_ylim(ep_ylow, ep_yhigh) + ax_ep.set_xlabel("时间 (s)") + ax_ep.set_ylabel("能量 / 功率") + ax_ep.set_title("系统能量与输入功率") + ax_ep.grid(True, alpha=0.3) + + line_ek, = ax_ep.plot([], [], "b-", lw=1.5, label="动能") + line_us, = ax_ep.plot([], [], "orange", lw=1.5, label="弹性势能") + line_et, = ax_ep.plot([], [], "r--", lw=1.5, label="总能量") + line_pw, = ax_ep.plot([], [], "g-", lw=1.5, alpha=0.7, label="输入功率 (dE/dt)") + if gravity_field: + line_ug, = ax_ep.plot([], [], "purple", lw=1.0, alpha=0.5, label="重力势能") + else: + line_ug = None + if gravity_interaction and n_atoms <= 200: + line_ugr, = ax_ep.plot([], [], "brown", lw=1.0, alpha=0.5, label="万有引力势能") + else: + line_ugr = None + + ax_ep.legend(loc="upper right", fontsize=9) + + plt.tight_layout(rect=[0, 0, 1, 0.95]) + + # ── 动画更新 ── + def update(frame): + # 波形 + for i in range(3): + wave_lines[i].set_data(atom_idx, [dx, dy, dz][i][frame]) + time_texts[i].set_text(f"t = {t[frame]:.2f} s | 帧 {frame+1}/{n_frames}") + + # 能量(累计到当前帧) + cur_t = t[:frame + 1] + line_ek.set_data(cur_t, ek[:frame + 1]) + line_us.set_data(cur_t, us[:frame + 1]) + line_et.set_data(cur_t, e_total[:frame + 1]) + line_pw.set_data(cur_t, power[:frame + 1]) + if line_ug: + line_ug.set_data(cur_t, ug[:frame + 1]) + if line_ugr: + line_ugr.set_data(cur_t, ugr[:frame + 1]) + + # 能量图 x 轴动态扩展 + ax_ep.set_xlim(t[0], max(t[frame] + max(t[-1] * 0.05, 1), t[-1])) + + artists = wave_lines + time_texts + [line_ek, line_us, line_et, line_pw] + if line_ug: artists.append(line_ug) + if line_ugr: artists.append(line_ugr) + return artists + + ani = FuncAnimation(fig, update, frames=n_frames, interval=50, blit=True) + + # 保存 GIF + gif_path = os.path.join(output_dir, "wave_animation.gif") + ani.save(gif_path, writer="pillow", fps=min(20, max(1, n_frames // 5))) + print(f"[plot_wave] GIF 已保存: {gif_path}") + + # 尝试保存 MP4 + try: + mp4_path = os.path.join(output_dir, "wave_animation.mp4") + ani.save(mp4_path, writer="ffmpeg", fps=min(20, max(1, n_frames // 5))) + print(f"[plot_wave] MP4 已保存: {mp4_path}") + except Exception: + pass + + plt.show() + return gif_path + + +if __name__ == "__main__": + script_dir = os.path.dirname(os.path.abspath(__file__)) + if len(sys.argv) > 1: + output_dir = os.path.abspath(sys.argv[1]) + else: + output_dir = compute.get_output_dir(script_dir) + plot_wave(output_dir)