求助！！！！调试程序时发现诡异现象！！！！ - Fortran Coder

[Fortran] 纯文本查看 复制代码

Subroutine checklimits

  Use param

  Implicit None

  ni1 = ni + 1

  niminx = -0.5*(lxx-llx) - 0.5*dx

  nimaxx = -0.5*dx

  niminz = -0.5*dz

  nimaxz = llz - 0.5*dz

  nominx = llx - 0.5*dx

  nomaxx = 0.5*(lxx+llx) - 0.5*dx

  nominz = -0.5*dz

  nomaxz = llz - 0.5*dz

  npminx = -0.5*dx

  npmaxx = llx - 0.5*dx

  npminz = -0.5*dz

  npmaxz = llz - 0.5*dz

  n2minx = -0.5*dx

  n2maxx = -0.5*dx + 0.5*(lxx-llx)

  n3minx = llx - 0.5*dx - 0.5*(lxx-llx)

  n3maxx = llx - 0.5*dx

  n5minxx = -0.5*dx + 0.5*(lxx-llx)

  n5maxx = llx - 0.5*dx - 0.5*(lxx-llx)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  Do k = 1, npmax

    xptemp(k) = 0.0

    zptemp(k) = 0.0

    uptemp(k) = 0.0

    wptemp(k) = 0.0

    ctemp(k) = 0.0

  End Do

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!check前，2,3,5三个区域各自的粒子数目分别为g2，g3，g5

  qq = 0

  no1 = no + 1 !!!!!此处修改过，9月9日

  Do k = no1, np

    If (xp(k)>=n2minx .And. xp(k)<=n2maxx) Then

      qq = qq + 1

    End If

  End Do

  g2 = qq

  qq = 0

  Do k = no1, np

    If (xp(k)>=n3minx .And. xp(k)<=n3maxx) Then

      qq = qq + 1

    End If

  End Do

  g3 = qq

  qq = 0

  Do k = no1, np

    If (xp(k)>=n5minxx .And. xp(k)<=n5maxx) Then

      qq = qq + 1

    End If

  End Do

  g5 = qq

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!更正出口粒子数no，更正流体粒子数np

  notemp = no + g3 - r

  nptemp = notemp + r + g5 + g3

  n21 = notemp + 1

  n22 = notemp + r

  n51 = notemp + r + 1

  n32 = notemp + r + g5 + g3

  n41 = ni + 1

  n42 = ni + g3

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!lv和lm的初始化

  Print *, 'LXX=', lxx

  Do i = 1, nptemp

    lv(i) = 0.0

    lm(i) = 0.0

  End Do

  Print *, 'LXX=', lxx

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!重新铺设2区域粒子,用临时转换标志temp表示!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  Do k = n21, n22

    xptemp(k) = xp(k-g3-r) + 0.5*(lxx-llx) !!!!!!!!!!!!!!!!!

    zptemp(k) = zp(k-g3-r)

    uptemp(k) = up(k-g3-r)

    wptemp(k) = wp(k-g3-r)

    ctemp(k) = c(k-g3-r)

  End Do

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!重新铺设4区域粒子!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  b = n41

  Do k = no1, np

    If (xp(k)>=n3minx .And. xp(k)<=n3maxx) Then

      lv(b) = k

      xpp(k) = xp(k) + 0.5*(lxx-llx)

      zpp(k) = zp(k)

      upp(k) = up(k)

      wpp(k) = wp(k)

      cpp(k) = c(k)

      b = b + 1

    End If

  End Do

  Do k = n41, n42

    xptemp(k) = xpp(lv(k))

    zptemp(k) = zpp(lv(k))

    uptemp(k) = upp(lv(k))

    wptemp(k) = wpp(lv(k))

    ctemp(k) = cpp(lv(k))

  End Do

!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!重新铺设53区域粒子!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  b = n51

  Do k = no1, np

    If (xp(k)>=n5minxx .And. xp(k)<=n3maxx) Then

      lm(b) = k

      xpp(k) = xp(k)

      zpp(k) = zp(k)

      upp(k) = up(k)

      wpp(k) = wp(k)

      cpp(k) = c(k)

      b = b + 1

    End If

  End Do

  Do k = n51, n32

    xptemp(k) = xpp(lm(k))

    zptemp(k) = zpp(lm(k))

    uptemp(k) = upp(lm(k))

    wptemp(k) = wpp(lm(k))

    ctemp(k) = cpp(lm(k))

  End Do

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  no = notemp

  np = nptemp

  Do k = ni1, np

    xp(k) = xptemp(k)

    zp(k) = zptemp(k)

    up(k) = uptemp(k)

    wp(k) = wptemp(k)

    c(k) = ctemp(k)

  End Do

End Subroutine checklimits

[Fortran] 纯文本查看 复制代码

Module param

  Implicit None

  Integer :: i, j, no1, np, kind_p2, lz, ncz, lx, ncx, lz2, nct, ii, jj, qq, no, g2, g3, g5, notemp, nptemp, r, nfile

  Integer :: n21, n22, n51, n32, n41, n42, ni, k, b, icell, kcell, nn, nnx, nnz, nbfd, nnn, nbfu, jj_start, time1

  Integer, Parameter :: npmax = 15000, kindp = 2, nct_max = 3000, nplink_max = 120, b1max = 100, b2max = 9940

  Integer :: lv(b1max), lm(b2max), iflag(npmax), ibox(nct_max, 2, nplink_max), nc(nct_max, kindp)

  Double Precision :: ax(npmax), az(npmax), ac(npmax), udot(npmax), wdot(npmax)

  Double Precision :: cdot(npmax), xp(npmax), zp(npmax), up(npmax)

  Double Precision :: wp(npmax), c(npmax), mass(npmax), xpp(b2max), zpp(b2max), upp(b2max), wpp(b2max)

  Double Precision :: xptemp(npmax), zptemp(npmax), uptemp(npmax), wptemp(npmax), ctemp(npmax)

  Double Precision :: cpp(b2max), xdot(npmax), zdot(npmax), xo(npmax), zo(npmax), uo(npmax), wo(npmax), co(npmax)

  Double Precision :: etap, d, drx, drz, rr2, fourh2, dux, duz, dcij, rr, hsml, tdwdr, eta, rho, niminx, nimaxx, niminz, nimaxz, nominx

  Double Precision :: nomaxx, nominz, nomaxz, npminx, npmaxx, npminz, npmaxz, n2minx, n2maxx, n3minx, n3maxx, n5minxx, n5maxx

  Double Precision :: xmin_ini, zmin_ini, xmax_ini, zmax_ini, one_over_2h, one_over_h, tw, factor, pi, g, time, dt, dt2, tmax

  Double Precision :: ncall11, ncall12, ncall13, h2

  Integer :: j1, j2, kind_p1, ini_kind_p2, lx2, ni1, n_start, n_end, kind_p

  Double Precision :: q

  Double Precision :: llx, llz, dx, dz, lxx, lzz, ddx, ddz

End Module param

Program sph

  Use param

  Implicit None

  ni = 4036

  ni1 = 4037

  llx = 2.5E-2

  llz = 2.5E-4

  lxx = 2.5E-2 + 4*2.5E-5

  lzz = 2.5E-4 + 4*2.5E-5

  dx = 2.5E-5

  dz = 2.5E-5

  no = 4056

  np = 14056

  r = 20

  k = 0

  Do i = 1, 14

    Do j = 1, 1004

      k = k + 1

      xp(nn) = (j-1)*dx - 0.5*(lxx-llx)

      zp(nn) = (i-1)*dz - 0.5*(lzz-llz)

      up(nn) = 0.0

      wp(nn) = 0.0

      c(nn) = 0.0

    End Do

  End Do

  Call checklimits

End Program sph