fix: Fortran engine use_marker/alpha header passthrough + case04 Moon orbit fix

- engines/fortran/main.f90:
  - 新增 use_marker 和 alpha_str 的读取、传递、写入 display.txt header
  - 新增 json_get_alpha 函数兼容 JSON 单值或 6 元数组
- engines/release/: 引擎校准缓存 & fortran.json 参数
- examples/case04/: 修正月球位置 (11.5→10.5) 和速度 (526→647),
  地月距缩至希尔半径以内 (0.5 < 1.0), 新增希尔半径验证说明
- examples/Readme.html: case04/case09 描述更新
This commit is contained in:
2026-06-19 08:05:05 +08:00
parent c8d12e7f45
commit 6c06c00f4b
8 changed files with 158 additions and 35 deletions
+70 -7
View File
@@ -21,6 +21,8 @@ program dynamics_f90
integer :: driving_force integer :: driving_force
double precision :: gravity_strength double precision :: gravity_strength
character(len=32) :: method character(len=32) :: method
integer :: use_marker
character(len=512) :: alpha_str
double precision :: t0, t1, elapsed, tw double precision :: t0, t1, elapsed, tw
! 原子数据 ! 原子数据
@@ -71,7 +73,8 @@ program dynamics_f90
! 读取 param.json ! 读取 param.json
call read_params(param_path, box_a, NT, DT, NSTEP, warmup_steps, method, G, B, & call read_params(param_path, box_a, NT, DT, NSTEP, warmup_steps, method, G, B, &
gravity_field, gravity_interaction, & 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 ! 读取 coord.txt
call read_coord(input_dir, n_atoms, atom_ids, masses, radii, pos_0, vel_0, fixed) 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, & traj_x, traj_y, traj_z, traj_vx, traj_vy, traj_vz, &
NT, DT, NSTEP, warmup_steps, method, G, B, & NT, DT, NSTEP, warmup_steps, method, G, B, &
n_bonds, gravity_field, elastic_force, damping_force, & 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) call cpu_time(t1)
elapsed = t1 - t0 elapsed = t1 - t0
@@ -228,15 +232,17 @@ contains
! ======================================================================== ! ========================================================================
subroutine read_params(path, box_a, NT, DT, NSTEP, warmup_steps, method, G, B, & subroutine read_params(path, box_a, NT, DT, NSTEP, warmup_steps, method, G, B, &
gravity_field, gravity_interaction, & 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 character(len=*), intent(in) :: path
double precision, intent(out) :: box_a, DT, G(3), B(3), gravity_strength double precision, intent(out) :: box_a, DT, G(3), B(3), gravity_strength
integer, intent(out) :: NT, NSTEP, warmup_steps integer, intent(out) :: NT, NSTEP, warmup_steps
integer, intent(out) :: gravity_field, gravity_interaction, elastic_force, damping_force, driving_force 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=8096) :: buf
character(len=256) :: line 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 box_a = 10.0d0; NT = 10000; DT = 0.001d0; NSTEP = 100; warmup_steps = 0
method = 'leapfrog' 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 gravity_field = 1; gravity_interaction = 0
elastic_force = 1; damping_force = 0; gravity_strength = 1.0d0 elastic_force = 1; damping_force = 0; gravity_strength = 1.0d0
driving_force = 0 driving_force = 0
use_marker = 0
alpha_str = '0.2'
open(newunit=u, file=trim(path), status='old', action='read', iostat=ios) open(newunit=u, file=trim(path), status='old', action='read', iostat=ios)
if (ios /= 0) then 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) damping_force = json_get_int(buf, 'damping_force', 0)
gravity_strength = json_get_double(buf, 'gravity_strength', 1.0d0) gravity_strength = json_get_double(buf, 'gravity_strength', 1.0d0)
driving_force = json_get_int(buf, 'driving_force', 0) 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 end subroutine read_params
! ======================================================================== ! ========================================================================
@@ -323,6 +333,55 @@ subroutine json_get_double3(buf, key, vals)
end do end do
end subroutine json_get_double3 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 中读取字符串值 ! 从 JSON 中读取字符串值
subroutine json_get_string(buf, key, val) subroutine json_get_string(buf, key, val)
character(len=*), intent(in) :: buf, key 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, & tx, ty, tz, tvx, tvy, tvz, &
NT, DT, NSTEP, warmup, method, G, B, & NT, DT, NSTEP, warmup, method, G, B, &
nb, gravity_field, elastic_force, damping_force, & nb, gravity_field, elastic_force, damping_force, &
driving_force, box_a, gravity_strength) driving_force, box_a, gravity_strength, &
character(len=*), intent(in) :: outdir, method 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) :: n_frames, nat, NT, NSTEP, warmup, nb
integer, intent(in) :: gravity_field, elastic_force, damping_force, driving_force integer, intent(in) :: gravity_field, elastic_force, damping_force, driving_force
integer, intent(in) :: use_marker
integer, intent(in) :: aid(nat) 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) :: 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) 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, '("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_MIN: ", g0)') -box_a; write(u, '(a)') trim(buf)
write(buf, '("Z_MAX: ", 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 ──────────────────────────────────────────────────────── ! ── frame data ────────────────────────────────────────────────────────
do f = 1, n_frames do f = 1, n_frames
+1 -1
View File
@@ -1 +1 @@
{"n_atoms": 40, "nt": 10000, "step_time": 0.00015027170181274412} {"n_atoms": 40, "nt": 10000, "step_time": 3.113259077072144e-05}
+1
View File
@@ -0,0 +1 @@
{"n_atoms": 40, "nt": 200000, "step_time": 0.0002908185601234436}
+51
View File
@@ -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
}
+6 -6
View File
@@ -210,7 +210,7 @@
<div class="case-card"> <div class="case-card">
<span class="num">04</span> <span class="num">04</span>
<h3>日地月系统(稳定)</h3> <h3>日地月系统(稳定)</h3>
<p>三体系统稳定椭圆轨道</p> <p>三体系统稳定轨道,经希尔半径检验</p>
<div class="meta"> <div class="meta">
<span class="tag tag-python">Python</span> <span class="tag tag-python">Python</span>
<span class="tag tag-gravity">万有引力</span> <span class="tag tag-gravity">万有引力</span>
@@ -263,7 +263,7 @@
<div class="case-card"> <div class="case-card">
<span class="num">09</span> <span class="num">09</span>
<h3>一维链纵波·Fortran 引擎</h3> <h3>一维链纵波·Fortran 引擎</h3>
<p>Fortran 引擎驱动的纵波传播测试</p> <p>Fortran 引擎驱动的纵波GPU 实例化渲染</p>
<div class="meta"> <div class="meta">
<span class="tag tag-fortran">Fortran</span> <span class="tag tag-fortran">Fortran</span>
<span class="tag tag-spring">弹簧</span> <span class="tag tag-spring">弹簧</span>
@@ -323,8 +323,8 @@
<li class="tag tag-python">Python 引擎</li> <li class="tag tag-python">Python 引擎</li>
<li class="tag tag-gravity">万有引力</li> <li class="tag tag-gravity">万有引力</li>
</ul> </ul>
<p>与 case03 相同的三体系统,但采用恰当的初始条件地球和月球维持稳定的椭圆轨道运动</p> <p>与 case03 相同的三体系统,但采用恰当的初始条件地球置于近日点($r=10$),$v_z=520$ 接近圆轨道;月球缩至希尔半径($r_H\approx1.0$)以内的 $r=0.5$,相对速度 $v_{\text{rel}}\approx127$,确保月球被地球稳定束缚</p>
<div class="highlight">✅ 与 case03 对比学习:初始条件对数值稳定性的影响</div> <div class="highlight">✅ 与 case03 对比学习:初始条件对数值稳定性的影响。地月距 0.5 在地球希尔半径以内,满足稳定性条件</div>
</div> </div>
<div class="detail-card"> <div class="detail-card">
@@ -377,8 +377,8 @@
<li class="tag tag-spring">弹簧键力</li> <li class="tag tag-spring">弹簧键力</li>
<li class="tag tag-drive">驱动力</li> <li class="tag tag-drive">驱动力</li>
</ul> </ul>
<p>40 原子沿 x 轴排列,Fortran 引擎驱动的纵波传播测试。验证 Fortran 引擎的输出兼容性。</p> <p>40 原子沿 x 轴排列,Fortran 引擎驱动的纵波传播测试。使用 <code>use_marker: 1</code>GPU 实例化 Marker 模式)加速渲染,验证 Fortran 引擎与其他引擎的输出兼容性。</p>
<div class="highlight">🔧 Fortran 引擎兼容性验证,T_total=200, NSTEP=100</div> <div class="highlight">🔧 Fortran 引擎兼容性验证,T_total=200, NSTEP=100。Marker 模式下 40 原子以 GPU 点精灵渲染,帧率高</div>
</div> </div>
<div class="detail-card"> <div class="detail-card">
+24 -16
View File
@@ -337,10 +337,14 @@ MathJax = {
<tr><th>天体</th><th>质量</th><th>位置 $(x,y,z)$</th><th>速度 $(v_x,v_y,v_z)$</th><th>固定约束</th></tr> <tr><th>天体</th><th>质量</th><th>位置 $(x,y,z)$</th><th>速度 $(v_x,v_y,v_z)$</th><th>固定约束</th></tr>
<tr><td>太阳</td><td>$27\,000$</td><td>$(0,0,0)$</td><td>$(0,0,0)$</td><td>全部固定</td></tr> <tr><td>太阳</td><td>$27\,000$</td><td>$(0,0,0)$</td><td>$(0,0,0)$</td><td>全部固定</td></tr>
<tr><td>地球</td><td>$81$</td><td>$(10,0,0)$</td><td>$(0,0,520)$</td><td></td></tr> <tr><td>地球</td><td>$81$</td><td>$(10,0,0)$</td><td>$(0,0,520)$</td><td></td></tr>
<tr><td>月球</td><td>$1$</td><td>$(11.5,0,0)$</td><td>$(0,0,526)$</td><td></td></tr> <tr><td>月球</td><td>$1$</td><td>$(10.5,0,0)$</td><td>$(0,0,647)$</td><td></td></tr>
</table> </table>
<p>地球在 $z$ 方向获得初速度 $v=520$,产生绕太阳的轨道运动;月球在地球基础上附加 $v=6$ 的绕地速度,形成绕地球的微小轨道。速度由圆形轨道公式 $v = \sqrt{G_{\text{eff}} M / r}$ 计算,其中 $G_{\text{eff}} = \text{gravity\_strength} = 100$。</p> <p>地球在 $z$ 方向获得初速度 $v=520$,产生绕太阳的轨道运动(接近圆轨道);月球在地球基础上附加 $v\approx127$ 的绕地速度,形成绕地球的轨道。速度由圆形轨道公式 $v = \sqrt{G_{\text{eff}} M / r}$ 计算,其中 $G_{\text{eff}} = \text{gravity\_strength} = 100$。</p>
<div class="highlight">
<strong>验证</strong>:地球速度 $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$ 之内,可确保轨道稳定。
</div>
<h2>五、使用方法</h2> <h2>五、使用方法</h2>
@@ -372,9 +376,10 @@ MathJax = {
<ul> <ul>
<li>太阳真实质量为月球 $27\,100\,000$ 倍,模拟中缩至 $27\,000$ 倍以保持数值稳定($1/1000$)</li> <li>太阳真实质量为月球 $27\,100\,000$ 倍,模拟中缩至 $27\,000$ 倍以保持数值稳定($1/1000$)</li>
<li>轨道半径未按真实比例缩放(真实日地距是地月距的 389 倍,模拟中6.7 倍</li> <li>轨道半径未按真实比例缩放(真实日地距是地月距的 389 倍,模拟中20 倍,受希尔半径约束,月球已尽可能放远</li>
<li>未考虑月球轨道倾角(真实地月轨道有约 5° 的倾角)</li> <li>未考虑月球轨道倾角(真实地月轨道有约 5° 的倾角)</li>
<li>leapfrog 算法为能量守恒的辛积分器,长期稳定,但步长过大时仍可能偏离真实轨道</li> <li>leapfrog 算法为能量守恒的辛积分器,长期稳定,但步长过大时仍可能偏离真实轨道</li>
<li>太阳固定不动可作为教学简化,但严格三体模拟中应让所有天体在质心系中自由运动</li>
</ul> </ul>
</div> </div>
@@ -424,18 +429,19 @@ MathJax = {
ctx.strokeStyle = '#c9d1d9'; ctx.strokeStyle = '#c9d1d9';
ctx.lineWidth = 1.5; ctx.lineWidth = 1.5;
ctx.stroke(); ctx.stroke();
// 箭头 // 箭头(左端:椭圆中心 cx+f
ctx.beginPath(); ctx.beginPath();
ctx.moveTo(cx - f + 6, aY - 4); ctx.moveTo(cx + f + 6, aY - 4);
ctx.lineTo(cx - f, aY); ctx.lineTo(cx + f, aY);
ctx.lineTo(cx - f + 6, aY + 4); ctx.lineTo(cx + f + 6, aY + 4);
ctx.strokeStyle = '#c9d1d9'; ctx.strokeStyle = '#c9d1d9';
ctx.lineWidth = 1.5; ctx.lineWidth = 1.5;
ctx.stroke(); ctx.stroke();
// 箭头(右端:右顶点 cx+f+A)
ctx.beginPath(); ctx.beginPath();
ctx.moveTo(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, aY);
ctx.lineTo(cx - f + A - 6, aY + 4); ctx.lineTo(cx + f + A - 6, aY + 4);
ctx.strokeStyle = '#c9d1d9'; ctx.strokeStyle = '#c9d1d9';
ctx.lineWidth = 1.5; ctx.lineWidth = 1.5;
ctx.stroke(); ctx.stroke();
@@ -453,17 +459,19 @@ MathJax = {
ctx.strokeStyle = '#d29922'; ctx.strokeStyle = '#d29922';
ctx.lineWidth = 1.5; ctx.lineWidth = 1.5;
ctx.stroke(); ctx.stroke();
// 箭头(左端:左焦点=cx 太阳)
ctx.beginPath(); ctx.beginPath();
ctx.moveTo(cx - f + 5, cY - 4); ctx.moveTo(cx + 5, cY - 4);
ctx.lineTo(cx - f, cY); ctx.lineTo(cx, cY);
ctx.lineTo(cx - f + 5, cY + 4); ctx.lineTo(cx + 5, cY + 4);
ctx.strokeStyle = '#d29922'; ctx.strokeStyle = '#d29922';
ctx.lineWidth = 1.5; ctx.lineWidth = 1.5;
ctx.stroke(); ctx.stroke();
// 箭头(右端:椭圆中心 cx+f)
ctx.beginPath(); ctx.beginPath();
ctx.moveTo(cx - 5, cY - 4); ctx.moveTo(cx + f - 5, cY - 4);
ctx.lineTo(cx, cY); ctx.lineTo(cx + f, cY);
ctx.lineTo(cx - 5, cY + 4); ctx.lineTo(cx + f - 5, cY + 4);
ctx.strokeStyle = '#d29922'; ctx.strokeStyle = '#d29922';
ctx.lineWidth = 1.5; ctx.lineWidth = 1.5;
ctx.stroke(); ctx.stroke();
+3 -3
View File
@@ -1,4 +1,4 @@
n mass radius x y z vx vy vz fix_x fix_y fix_z 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 1 27000 1.0 0 0 0 0 0 0 1 1 1
2 81 0.5 10 0 0 0 0 520 0 0 0 2 81 0.2 10 0 0 0 0 520 0 0 0
3 1 0.3 11.5 0 0 0 0 526 0 0 0 3 1 0.1 10.5 0 0 0 0 647 0 0 0
+1 -1
View File
@@ -74,7 +74,7 @@ T_total: 10.0
NSTEP: 2 NSTEP: 2
# ── 时间步长 ────────────────────────────────── # ── 时间步长 ──────────────────────────────────
DT: 0.001 # 时间步长 (s) DT: 0.0001 # 时间步长 (s)
# 抽帧范围:只保存 [sample_start, sample_end) 区间内的帧 # 抽帧范围:只保存 [sample_start, sample_end) 区间内的帧
sample_start: null # null 表示从头开始(帧索引从 0 起) sample_start: null # null 表示从头开始(帧索引从 0 起)