diff --git a/engines/fortran/main.f90 b/engines/fortran/main.f90 index 2e90c05..06a8afc 100644 --- a/engines/fortran/main.f90 +++ b/engines/fortran/main.f90 @@ -21,6 +21,8 @@ program dynamics_f90 integer :: driving_force double precision :: gravity_strength character(len=32) :: method + integer :: use_marker + character(len=512) :: alpha_str double precision :: t0, t1, elapsed, tw ! 原子数据 @@ -71,7 +73,8 @@ program dynamics_f90 ! 读取 param.json call read_params(param_path, box_a, NT, DT, NSTEP, warmup_steps, method, G, B, & gravity_field, gravity_interaction, & - elastic_force, damping_force, gravity_strength, driving_force) + elastic_force, damping_force, gravity_strength, driving_force, & + use_marker, alpha_str) ! 读取 coord.txt call read_coord(input_dir, n_atoms, atom_ids, masses, radii, pos_0, vel_0, fixed) @@ -201,7 +204,8 @@ program dynamics_f90 traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz, & NT, DT, NSTEP, warmup_steps, method, G, B, & n_bonds, gravity_field, elastic_force, damping_force, & - driving_force, box_a, gravity_strength) + driving_force, box_a, gravity_strength, & + use_marker, alpha_str) call cpu_time(t1) elapsed = t1 - t0 @@ -228,15 +232,17 @@ contains ! ======================================================================== subroutine read_params(path, box_a, NT, DT, NSTEP, warmup_steps, method, G, B, & gravity_field, gravity_interaction, & - elastic_force, damping_force, gravity_strength, driving_force) + elastic_force, damping_force, gravity_strength, driving_force, & + use_marker, alpha_str) character(len=*), intent(in) :: path double precision, intent(out) :: box_a, DT, G(3), B(3), gravity_strength integer, intent(out) :: NT, NSTEP, warmup_steps integer, intent(out) :: gravity_field, gravity_interaction, elastic_force, damping_force, driving_force - character(len=*), intent(out) :: method + integer, intent(out) :: use_marker + character(len=*), intent(out) :: method, alpha_str character(len=8096) :: buf character(len=256) :: line - integer :: u, ios + integer :: u, ios, i box_a = 10.0d0; NT = 10000; DT = 0.001d0; NSTEP = 100; warmup_steps = 0 method = 'leapfrog' @@ -244,6 +250,8 @@ subroutine read_params(path, box_a, NT, DT, NSTEP, warmup_steps, method, G, B, & gravity_field = 1; gravity_interaction = 0 elastic_force = 1; damping_force = 0; gravity_strength = 1.0d0 driving_force = 0 + use_marker = 0 + alpha_str = '0.2' open(newunit=u, file=trim(path), status='old', action='read', iostat=ios) if (ios /= 0) then @@ -274,6 +282,8 @@ subroutine read_params(path, box_a, NT, DT, NSTEP, warmup_steps, method, G, B, & damping_force = json_get_int(buf, 'damping_force', 0) gravity_strength = json_get_double(buf, 'gravity_strength', 1.0d0) driving_force = json_get_int(buf, 'driving_force', 0) + use_marker = json_get_int(buf, 'use_marker', 0) + call json_get_alpha(buf, 'alpha', alpha_str) end subroutine read_params ! ======================================================================== @@ -323,6 +333,55 @@ subroutine json_get_double3(buf, key, vals) end do end subroutine json_get_double3 +! 读取 alpha：兼容单值 (0.2) 或 6 元数组 [0.0,0.0,0.0,0.0,0.0,0.0] +! 输出为逗号分隔字符串，与 draw.py 的解析格式一致 +subroutine json_get_alpha(buf, key, val) + character(len=*), intent(in) :: buf, key + character(len=*), intent(out) :: val + double precision :: d6(6) + integer :: p, i, ios + character(len=32) :: tmp + val = '0.2' + p = index(buf, '"' // key // '"') + if (p == 0) return + p = index(buf(p:), ':') + p + ! skip whitespace + do while (p <= len(buf) .and. (buf(p:p) == ':' .or. buf(p:p) == ' ' & + .or. buf(p:p) == char(9) .or. buf(p:p) == char(10))) + p = p + 1 + end do + if (buf(p:p) == '[') then + ! 6 元数组 + d6 = 0.2d0 + p = p + 1 + do i = 1, 6 + do while (p <= len(buf) .and. (buf(p:p) == ' ' .or. buf(p:p) == ',' & + .or. buf(p:p) == char(9) .or. buf(p:p) == char(10))) + p = p + 1 + end do + if (p > len(buf) .or. buf(p:p) == ']') exit + read(buf(p:), *, iostat=ios) d6(i) + if (ios /= 0) exit + do while (p <= len(buf) .and. buf(p:p) /= ',' .and. buf(p:p) /= ']') + p = p + 1 + end do + end do + val = '' + do i = 1, 6 + write(tmp, '(f0.6)') d6(i) + if (i > 1) val = trim(val) // ',' + val = trim(val) // trim(tmp) + end do + else + ! 单值标量 + read(buf(p:), *, iostat=ios) d6(1) + if (ios == 0) then + write(tmp, '(f0.6)') d6(1) + val = trim(tmp) + end if + end if +end subroutine json_get_alpha + ! 从 JSON 中读取字符串值 subroutine json_get_string(buf, key, val) character(len=*), intent(in) :: buf, key @@ -999,10 +1058,12 @@ subroutine write_display_txt(outdir, n_frames, nat, aid, & tx, ty, tz, tvx, tvy, tvz, & NT, DT, NSTEP, warmup, method, G, B, & nb, gravity_field, elastic_force, damping_force, & - driving_force, box_a, gravity_strength) - character(len=*), intent(in) :: outdir, method + driving_force, box_a, gravity_strength, & + use_marker, alpha_str) + character(len=*), intent(in) :: outdir, method, alpha_str integer, intent(in) :: n_frames, nat, NT, NSTEP, warmup, nb integer, intent(in) :: gravity_field, elastic_force, damping_force, driving_force + integer, intent(in) :: use_marker integer, intent(in) :: aid(nat) double precision, intent(in) :: tx(n_frames, nat), ty(n_frames, nat), tz(n_frames, nat) double precision, intent(in) :: tvx(n_frames, nat), tvy(n_frames, nat), tvz(n_frames, nat) @@ -1051,6 +1112,8 @@ subroutine write_display_txt(outdir, n_frames, nat, aid, & write(buf, '("Y_MAX: ", g0)') box_a; write(u, '(a)') trim(buf) write(buf, '("Z_MIN: ", g0)') -box_a; write(u, '(a)') trim(buf) write(buf, '("Z_MAX: ", g0)') box_a; write(u, '(a)') trim(buf) + write(u, '("use_marker: ", i0)') use_marker + write(u, '("alpha: ", a)') trim(alpha_str) ! ── frame data ──────────────────────────────────────────────────────── do f = 1, n_frames diff --git a/engines/release/_calib_c.json b/engines/release/_calib_c.json index d19b77c..c879a09 100644 --- a/engines/release/_calib_c.json +++ b/engines/release/_calib_c.json @@ -1 +1 @@ -{"n_atoms": 40, "nt": 10000, "step_time": 0.00015027170181274412} \ No newline at end of file +{"n_atoms": 40, "nt": 10000, "step_time": 3.113259077072144e-05} \ No newline at end of file diff --git a/engines/release/_calib_fortran.json b/engines/release/_calib_fortran.json new file mode 100644 index 0000000..3ec89e5 --- /dev/null +++ b/engines/release/_calib_fortran.json @@ -0,0 +1 @@ +{"n_atoms": 40, "nt": 200000, "step_time": 0.0002908185601234436} \ No newline at end of file diff --git a/engines/release/fortran.json b/engines/release/fortran.json new file mode 100644 index 0000000..a25e603 --- /dev/null +++ b/engines/release/fortran.json @@ -0,0 +1,51 @@ +{ + "box_a": 300.0, + "NT": 200000, + "DT": 0.001, + "NSTEP": 100, + "warmup_steps": 0, + "method": "leapfrog", + "G": [ + 0.0, + 0.0, + 0.0 + ], + "B": [ + 0.005, + 0.0, + 0.005 + ], + "gravity_field": 0, + "gravity_interaction": 0, + "elastic_force": 1, + "damping_force": 0, + "gravity_strength": 1.0, + "driving_force": 1, + "save_trajectory": 0, + "alpha": [ + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0 + ], + "ball_radius": 0.5, + "ball_color": [ + 0.2, + 0.6, + 0.9 + ], + "box_color": [ + 0.8, + 0.8, + 0.85 + ], + "use_marker": 1, + "camera_distance": 60.0, + "camera_elevation": 0.0, + "camera_azimuth": 0.0, + "camera_center_x": 30.0, + "camera_center_y": 0.0, + "camera_center_z": 0.0 +} \ No newline at end of file diff --git a/examples/Readme.html b/examples/Readme.html index 5aa9d2b..8443a78 100644 --- a/examples/Readme.html +++ b/examples/Readme.html @@ -210,7 +210,7 @@
04

日地月系统（稳定）

-

三体系统稳定椭圆轨道

+

三体系统稳定轨道，经希尔半径检验

Python 万有引力 @@ -263,7 +263,7 @@
09

一维链纵波·Fortran 引擎

-

Fortran 引擎驱动的纵波传播测试

+

Fortran 引擎驱动的纵波，GPU 实例化渲染

Fortran 弹簧 @@ -323,8 +323,8 @@
  • Python 引擎
  • 万有引力
    • -

    与 case03 相同的三体系统，但采用恰当的初始条件，地球和月球维持稳定的椭圆轨道运动。

    -
    ✅ 与 case03 对比学习：初始条件对数值稳定性的影响
    +

    与 case03 相同的三体系统，但采用恰当的初始条件：地球置于近日点（$r=10$），$v_z=520$ 接近圆轨道；月球缩至希尔半径（$r_H\approx1.0$）以内的 $r=0.5$，相对速度 $v_{\text{rel}}\approx127$，确保月球被地球稳定束缚。

    +
    ✅ 与 case03 对比学习：初始条件对数值稳定性的影响。地月距 0.5 在地球希尔半径以内，满足稳定性条件
    @@ -377,8 +377,8 @@
  • 弹簧键力
  • 驱动力
    • -

    40 原子沿 x 轴排列，Fortran 引擎驱动的纵波传播测试。验证 Fortran 引擎的输出兼容性。

    -
    🔧 Fortran 引擎兼容性验证，T_total=200, NSTEP=100
    +

    40 原子沿 x 轴排列，Fortran 引擎驱动的纵波传播测试。使用 use_marker: 1（GPU 实例化 Marker 模式）加速渲染，验证 Fortran 引擎与其他引擎的输出兼容性。

    +
    🔧 Fortran 引擎兼容性验证，T_total=200, NSTEP=100。Marker 模式下 40 原子以 GPU 点精灵渲染，帧率高
    diff --git a/examples/case04/Readme.html b/examples/case04/Readme.html index 4c2e8cf..2337a83 100644 --- a/examples/case04/Readme.html +++ b/examples/case04/Readme.html @@ -337,10 +337,14 @@ MathJax = { 天体质量位置 $(x,y,z)$速度 $(v_x,v_y,v_z)$固定约束 太阳$27\,000$$(0,0,0)$$(0,0,0)$全部固定 地球$81$$(10,0,0)$$(0,0,520)$无 - 月球$1$$(11.5,0,0)$$(0,0,526)$无 + 月球$1$$(10.5,0,0)$$(0,0,647)$无 -

    地球在 $z$ 方向获得初速度 $v=520$，产生绕太阳的轨道运动；月球在地球基础上附加 $v=6$ 的绕地速度，形成绕地球的微小轨道。速度由圆形轨道公式 $v = \sqrt{G_{\text{eff}} M / r}$ 计算，其中 $G_{\text{eff}} = \text{gravity\_strength} = 100$。

    +

    地球在 $z$ 方向获得初速度 $v=520$，产生绕太阳的轨道运动（接近圆轨道）；月球在地球基础上附加 $v\approx127$ 的绕地速度，形成绕地球的轨道。速度由圆形轨道公式 $v = \sqrt{G_{\text{eff}} M / r}$ 计算，其中 $G_{\text{eff}} = \text{gravity\_strength} = 100$。

    + +
    + 验证：地球速度 $v_\oplus = \sqrt{100 \times 27\,000 / 10} \approx 519.6$，设 520 正确。月球相对速度 $v_{\text{rel}} = \sqrt{100 \times 81 / 0.5} \approx 127.3$，设 127 正确。地月距离 0.5 在地球希尔半径 $r_H \approx 10 \times (81 / 81\,000)^{1/3} \approx 1.0$ 之内，可确保轨道稳定。 +

    五、使用方法

    @@ -372,9 +376,10 @@ MathJax = {
    • 太阳真实质量为月球 $27\,100\,000$ 倍，模拟中缩至 $27\,000$ 倍以保持数值稳定（$1/1000$）
      • -
    • 轨道半径未按真实比例缩放（真实日地距是地月距的 389 倍，模拟中仅为 6.7 倍）
      • +
    • 轨道半径未按真实比例缩放（真实日地距是地月距的 389 倍，模拟中约为 20 倍，受希尔半径约束，月球已尽可能放远）
    • 未考虑月球轨道倾角（真实地月轨道有约 5° 的倾角）
    • leapfrog 算法为能量守恒的辛积分器，长期稳定，但步长过大时仍可能偏离真实轨道
      • +
    • 太阳固定不动可作为教学简化，但严格三体模拟中应让所有天体在质心系中自由运动
    @@ -424,18 +429,19 @@ MathJax = { ctx.strokeStyle = '#c9d1d9'; ctx.lineWidth = 1.5; ctx.stroke(); - // 箭头 + // 箭头（左端：椭圆中心 cx+f） ctx.beginPath(); - ctx.moveTo(cx - f + 6, aY - 4); - ctx.lineTo(cx - f, aY); - ctx.lineTo(cx - f + 6, aY + 4); + ctx.moveTo(cx + f + 6, aY - 4); + ctx.lineTo(cx + f, aY); + ctx.lineTo(cx + f + 6, aY + 4); ctx.strokeStyle = '#c9d1d9'; ctx.lineWidth = 1.5; ctx.stroke(); + // 箭头（右端：右顶点 cx+f+A） ctx.beginPath(); - ctx.moveTo(cx - f + A - 6, aY - 4); - ctx.lineTo(cx - f + A, aY); - ctx.lineTo(cx - f + A - 6, aY + 4); + ctx.moveTo(cx + f + A - 6, aY - 4); + ctx.lineTo(cx + f + A, aY); + ctx.lineTo(cx + f + A - 6, aY + 4); ctx.strokeStyle = '#c9d1d9'; ctx.lineWidth = 1.5; ctx.stroke(); @@ -453,17 +459,19 @@ MathJax = { ctx.strokeStyle = '#d29922'; ctx.lineWidth = 1.5; ctx.stroke(); + // 箭头（左端：左焦点=cx 太阳） ctx.beginPath(); - ctx.moveTo(cx - f + 5, cY - 4); - ctx.lineTo(cx - f, cY); - ctx.lineTo(cx - f + 5, cY + 4); + ctx.moveTo(cx + 5, cY - 4); + ctx.lineTo(cx, cY); + ctx.lineTo(cx + 5, cY + 4); ctx.strokeStyle = '#d29922'; ctx.lineWidth = 1.5; ctx.stroke(); + // 箭头（右端：椭圆中心 cx+f） ctx.beginPath(); - ctx.moveTo(cx - 5, cY - 4); - ctx.lineTo(cx, cY); - ctx.lineTo(cx - 5, cY + 4); + ctx.moveTo(cx + f - 5, cY - 4); + ctx.lineTo(cx + f, cY); + ctx.lineTo(cx + f - 5, cY + 4); ctx.strokeStyle = '#d29922'; ctx.lineWidth = 1.5; ctx.stroke(); diff --git a/examples/case04/input/coord.txt b/examples/case04/input/coord.txt index 036ade1..25f439d 100644 --- a/examples/case04/input/coord.txt +++ b/examples/case04/input/coord.txt @@ -1,4 +1,4 @@ -n mass radius x y z vx vy vz fix_x fix_y fix_z -1 27000 2.0 0 0 0 0 0 0 1 1 1 -2 81 0.5 10 0 0 0 0 520 0 0 0 -3 1 0.3 11.5 0 0 0 0 526 0 0 0 \ No newline at end of file +n mass radius x y z vx vy vz fix_x fix_y fix_z +1 27000 1.0 0 0 0 0 0 0 1 1 1 +2 81 0.2 10 0 0 0 0 520 0 0 0 +3 1 0.1 10.5 0 0 0 0 647 0 0 0 \ No newline at end of file diff --git a/examples/case04/input/input.txt b/examples/case04/input/input.txt index 1bf6fbd..61432a6 100644 --- a/examples/case04/input/input.txt +++ b/examples/case04/input/input.txt @@ -74,7 +74,7 @@ T_total: 10.0 NSTEP: 2 # ── 时间步长 ────────────────────────────────── -DT: 0.001 # 时间步长 (s) +DT: 0.0001 # 时间步长 (s) # 抽帧范围：只保存 [sample_start, sample_end) 区间内的帧 sample_start: null # null 表示从头开始（帧索引从 0 起）