以下是我本周末为一些介绍性研究工作而编写的一维水电代码(使用哈密顿方法和斯特朗 split 来演化变量 p 和 q)
do
if(num==1) then
p2 = p(i) - (dt/2.)*q(i)/abs(q(i)) ! half step in P
q(i) = q(i) + dt*p2 ! full step in Q
p(i) = p2 - (dt/2.)*q(i)/abs(q(i)) ! half step in P
num=2
elseif(num==2) then
q2 = q(i) + (dt/2.)*p(i) ! half step in Q
p(i) = p(i) - dt*q2/abs(q2) ! full step in P
q(i) = q(i) + (dt/2.)*p(i) ! half step in Q
num=1
endif
t = t+dt
if(t >= tend) exit
enddo
有没有比我这里有的更有效的方法来交替使用这两种算法(这是减少虚假数据所必需的)?如果重要的话,p 和 q 各有大约 100,000 个单元格(代码是并行化的)。
编辑:我添加了代码的 do
-loop 部分,而不仅仅是 if-elseif
部分。在 endif
之后还有一个写入文件的部分,但我认为这对于潜在的优化不是必需的。
最佳答案
我会重写代码以完全删除 if/then/else
:
integer :: num_steps, k
logical :: one_more
num_steps = tend/dt
one_more = (mod(num_steps,2) /= 0)
do k = 1,num_steps/2
p2 = p(i) - (dt/2.)*q(i)/abs(q(i)) ! half step in P
q(i) = q(i) + dt*p2 ! full step in Q
p(i) = p2 - (dt/2.)*q(i)/abs(q(i)) ! half step in P
! output
q2 = q(i) + (dt/2.)*p(i) ! half step in Q
p(i) = p(i) - dt*q2/abs(q2) ! full step in P
q(i) = q(i) + (dt/2.)*p(i) ! half step in Q
! output
enddo
if (one_more) then
p2 = p(i) - (dt/2.)*q(i)/abs(q(i)) ! half step in P
q(i) = q(i) + dt*p2 ! full step in Q
p(i) = p2 - (dt/2.)*q(i)/abs(q(i)) ! half step in P
! output
endif
如果您需要当前时间进行输出操作,您仍然可以在循环中的每个步骤之后使用 t = t+dt
语句。
关于algorithm - 交替拆分运算符方法,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/13674754/