Optimized: Expand loops in dynamics/adv/ffsl_mod.F90loop

[shepherdmrf]

Written in #24

[shepherdmrf]

[shepherdmrf]

I rewrite the if-else part , please help me check my code. Thank you! @BeverlyCrl

[BeverlyCrl]

咱觉得这个if-else part有点奇怪了。还是分支的问题，如下图，你是这样处理的。 先判断i和i+1是否都置零，如果没有都置零进入最外层else。 进入这个else有三种情况：i置零，i+1置零，两者都不置零 考虑到此时两者同时置零不可能发生，最前面的if和else if是🆗的。然后就进入了第二张图里面的里层else部分，请看下面。

这里让咱先贴一下最原始的代码，这里的if-else结构可以看出有这么3个分支：数字极小时置零，如果没有置零且为正数进elif，剩下的情况进else。 也就是说 绝对值够小进if 不够小时正的进elif，负的进else

然后让我们看看这里你修改后的部分。此时能够进入else部分的只有一种情况，i和i+1都没有被置零。然后在这个部分你处理了i和i+1的原代码中elif和else部分。但是存在这样一种情况，当只有一方置零，比如i的置零了。此时i+1应该单独于i考虑进入原代码中的elif或者else部分。但是在你的代码中因为已经进入if (abs(cflx%d(i,j,k)) < 1.0e-16_r8) then这个分支中，不会在进入else中处理了。这个时候的i+1是没有被处理的。 当i+1置零，i不置零时同理。

这里面来说 你同时想要测试两个策略在这里的优化，会有一些confused其实。 可以拆开来尝试，其实我最开始是想让你试试只有循环展开没有改变任何ifelse的时候 这个代码会不会有效果。如果有效果 在此基础上再进行ifelse的重构 再比较可能会让结果清晰很多。纯粹的循环展开其实很简单也很笨拙，如下面的代码（当然只展开2的话最后完全可以写成if）：

        do i = mesh%half_ids, mesh%half_ide - 1, 2
            ci = int(cflx%d(i,j,k))
            cf = cflx%d(i,j,k) - ci
            if (abs(cflx%d(i,j,k)) < 1.0e-16_r8) then
              mfx%d(i,j,k) = 0
            else if (cflx%d(i,j,k) > 0) then
              iu = i - ci
              call ppm(mx%d(iu-2,j,k), mx%d(iu-1,j,k), mx%d(iu,j,k), mx%d(iu+1,j,k), mx%d(iu+2,j,k), ml, dm, m6)
              s1 = 1 - cf
              s2 = 1
              ds1 = s2    - s1
              ds2 = s2**2 - s1**2
              ds3 = s2**3 - s1**3
              mfx%d(i,j,k) =  u%d(i,j,k) * (sum(mx%d(iu+1:i,j,k)) + ml * ds1 + 0.5_r8 * dm * ds2 + m6 * (ds2 / 2.0_r8 - ds3 / 3.0_r8)) / cflx%d(i,j,k)
            else
              iu = i - ci + 1
              call ppm(mx%d(iu-2,j,k), mx%d(iu-1,j,k), mx%d(iu,j,k), mx%d(iu+1,j,k), mx%d(iu+2,j,k), ml, dm, m6)
              s1 = 0
              s2 = -cf
              ds1 = s2    - s1
              ds2 = s2**2 - s1**2
              ds3 = s2**3 - s1**3
              mfx%d(i,j,k) = -u%d(i,j,k) * (sum(mx%d(i+1:iu-1,j,k)) + ml * ds1 + 0.5_r8 * dm * ds2 + m6 * (ds2 / 2.0_r8 - ds3 / 3.0_r8)) / cflx%d(i,j,k)
            end if

            ci = int(cflx%d(i + 1,j,k))
            cf = cflx%d(i + 1,j,k) - ci
            if (abs(cflx%d(i,j,k)) < 1.0e-16_r8) then
              mfx%d(i + 1,j,k) = 0
            else if (cflx%d(i + 1,j,k) > 0) then
              iu = i + 1 - ci
              call ppm(mx%d(iu-2,j,k), mx%d(iu-1,j,k), mx%d(iu,j,k), mx%d(iu+1,j,k), mx%d(iu+2,j,k), ml, dm, m6)
              s1 = 1 - cf
              s2 = 1
              ds1 = s2    - s1
              ds2 = s2**2 - s1**2
              ds3 = s2**3 - s1**3
              mfx%d(i + 1,j,k) =  u%d(i + 1,j,k) * (sum(mx%d(iu+1:i + 1,j,k)) + ml * ds1 + 0.5_r8 * dm * ds2 + m6 * (ds2 / 2.0_r8 - ds3 / 3.0_r8)) / cflx%d(i + 1,j,k)
            else
              iu = i + 1 - ci + 1
              call ppm(mx%d(iu-2,j,k), mx%d(iu-1,j,k), mx%d(iu,j,k), mx%d(iu+1,j,k), mx%d(iu+2,j,k), ml, dm, m6)
              s1 = 0
              s2 = -cf
              ds1 = s2    - s1
              ds2 = s2**2 - s1**2
              ds3 = s2**3 - s1**3
              mfx%d(i + 1,j,k) = -u%d(i,j,k) * (sum(mx%d(i+2:iu-1,j,k)) + ml * ds1 + 0.5_r8 * dm * ds2 + m6 * (ds2 / 2.0_r8 - ds3 / 3.0_r8)) / cflx%d(i + 1,j,k)
            end if

        end do

        do i = mesh%half_ide - mod(mesh%half_ide - mesh%half_ids, 2), mesh%half_ide
            ci = int(cflx%d(i,j,k))
            cf = cflx%d(i,j,k) - ci
            if (abs(cflx%d(i,j,k)) < 1.0e-16_r8) then
              mfx%d(i,j,k) = 0
            else if (cflx%d(i,j,k) > 0) then
              iu = i - ci
              call ppm(mx%d(iu-2,j,k), mx%d(iu-1,j,k), mx%d(iu,j,k), mx%d(iu+1,j,k), mx%d(iu+2,j,k), ml, dm, m6)
              s1 = 1 - cf
              s2 = 1
              ds1 = s2    - s1
              ds2 = s2**2 - s1**2
              ds3 = s2**3 - s1**3
              mfx%d(i,j,k) =  u%d(i,j,k) * (sum(mx%d(iu+1:i,j,k)) + ml * ds1 + 0.5_r8 * dm * ds2 + m6 * (ds2 / 2.0_r8 - ds3 / 3.0_r8)) / cflx%d(i,j,k)
            else
              iu = i - ci + 1
              call ppm(mx%d(iu-2,j,k), mx%d(iu-1,j,k), mx%d(iu,j,k), mx%d(iu+1,j,k), mx%d(iu+2,j,k), ml, dm, m6)
              s1 = 0
              s2 = -cf
              ds1 = s2    - s1
              ds2 = s2**2 - s1**2
              ds3 = s2**3 - s1**3
              mfx%d(i,j,k) = -u%d(i,j,k) * (sum(mx%d(i+1:iu-1,j,k)) + ml * ds1 + 0.5_r8 * dm * ds2 + m6 * (ds2 / 2.0_r8 - ds3 / 3.0_r8)) / cflx%d(i,j,k)
            end if
        end do

[shepherdmrf]

Thank you! I will do all my power to solve this