程序计算时，参数为什么老是为0呢？

cainiao07 · 发表于 2014-3-21 11:40:35

本人目前在做岩石损伤本构的二次开发，使用fotran编译。现在有一个问题比较恼火，就是
参数ran的随机数是可以得到的，但是通过 et_r=cm(5)*(-1.0*log(ran1))**(1.0/wb1)计算后就
变为0了

我尝试去改变计算式，发现不行，怀疑是随机程序的问题，使用fotran帮助用ranom_number()
计算的hsv(5)还是0.大神们帮我看看是什么原因，小弟万分感谢！

下面是监控的截图，图中两个0,都是上述公式计算的。但是ran的值是有的

c    随机产生切线模量,下限等于0以避免计算发散，切线模量的平均值为输入的cm(5)
c    第一次计算步抽样，然后存储到hisv(5)中。
   call seed(1995)
   if(tt.LE.dt1) then
   do 101 i=1,cm(11)
      CALL RANDOM(ran1) !取得随机的随机数，否则每次进入系统（重新计算），抽样值都相同
101    enddo
   endif
c    CALL RANDOM_SEED()
c    CALL RANDOM_NUMBER(HARVEST =RAN)
if(tt.LE.dt1) then
   et_r=cm(5)*(-1.0*log(ran1))**(1.0/wb1)  !hisv(5)即为服从Weibull分布
   hsv(5)=et_r
endif

楚香饭 · 发表于 2014-3-21 19:42:24

黑色的ran值有吗？e-35 就跟0一样。你确定是你想要的吗？试试换一个种子呢？

如果你确定ran值是有的，那就不必怀疑随机数了。看看 wb1,cm(5) 的值，及其类型吧。

cainiao07 · 发表于 2014-3-21 20:47:53

chuxf 发表于 2014-3-21 19:42
黑色的ran值有吗？e-35 就跟0一样。你确定是你想要的吗？试试换一个种子呢？

如果你确定ran值是有的，那就 ...

这个是源程序，其中产生随机数是看的网上大家推荐的那种，另一种被盖住的是以前老版本提供的一段程序。我在做二次开发，调试了两周，问题就出在这，麻烦帮我看一下，不甚感激。

程序中是一段weibul分布的随机抽样，问题就出在这里了l

c    call seed(1995)
   if(tt.LE.dt1) then
c    do 101 i=1,cm(11)
c       CALL RANDOM(ran1) !取得随机的随机数，否则每次进入系统（重新计算），抽样值都相同
c101    enddo
   CALL RANDOM_SEED()
   CALL RANDOM_NUMBER(HARVEST =RAN1)
   et_r=10000*(-1.0*log(ran1))**(1.0/wb1)
   !hisv(5)即为服从Weibull分布

   hsv(5)=et_r
c    write(*,*) hsv(5)
   endif

楚香饭 · 发表于 2014-3-21 21:31:54

1.源程序不完整，我无法帮助你调试。
2.我没做过二次开发，无法帮助你调试。
3.random_seed 只许调用一次，请放置在程序一开始（不能在循环内）
4.请跟踪调试产生的随机数，是多少？是否满足你的期望？
5.如随机数产生没有问题，后续计算出现0，则不是随机数的问题。请跟踪 wb1 和 cm(5)

cainiao07 · 发表于 2014-3-22 23:05:26

chuxf 发表于 2014-3-21 21:31
1.源程序不完整，我无法帮助你调试。
2.我没做过二次开发，无法帮助你调试。
3.random_seed 只许调用一次， ...

大哥你好，我把程序抠出来单独编译。为了与二次开发相吻合，取了相近的数值。
计算就发现WX_T,ET1没有数值,我实验过很多次，怀疑原因为：当前面有一段随机数导致的。

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

!****************************************************************************
          subroutine init_random_seed()
       integer :: i,n,clock
       integer,dimension(:),allocatable :: seed
       call random_seed(size=n)
        allocate(seed(n))
        call system_clock(count=clock)
        seed=clock+37*(/(i-1,i=1,n)/)
        call random_seed(PUT=seed) 
        deallocate(seed)
      end subroutine

    program Console2
    implicit none
    REAL WB_X,hsv1,hsv2,rd     
        REAL ::ET1
        real ::ran         

       call init_random_seed()   
       CALL RANDOM_NUMBER(ran)

          WB_X=10000*(-1.0*(ran))**(1.0/5)  !WB_X即为服从Weibull分布 的一个抽样
          ET1=10000
          ET1=ET1*(-1.0*(ran))**(1.0/5)
          hsv1=ET1                             !hisv(5)即为服从Weibull分布
          hsv2=WB_X 
                                  !初始损伤服从webull分布
       print(*,*),ET1,WB_X,RAN,HSV1,hsv2 
       stop 
    end program Console2

楚香饭 · 发表于 2014-3-23 12:29:49

WB_X=10000*(-1.0*(ran))**(1.0/5)

我建议你写为

WB_X=-10000*((ran))**(1.0/5)

虽然数学上负数允许开奇次方，但是由于 0.1/5 本身存在是浮点数，可能有误差，编译器无法准确判断一个数是否为奇数。

因此建议把负号，提到开方前面去。

后面的 ET1=ET1*(-1.0*(ran))**(1.0/5) 同理

cainiao07 · 发表于 2014-3-23 21:02:38

chuxf 发表于 2014-3-23 12:29
WB_X=10000*(-1.0*(ran))**(1.0/5)

我建议你写为

感谢您的不懈回复，虽然没有调试通，但是大致分析出了问题所在：
1、负数不能开奇次方
2、LOG函数在LS-DYNA源程序二次开发中属于外部函数调用，估计是和主程序冲突吧
3、我在程序中加入了n=n+1进行监视，发现整个本构程序在计算时候只进行了一次调用，说明程序在中间没有执行，具体原因不详，有待分析。但是至少发现了这个计算无值的根本原因。

楚香饭 · 发表于 2014-3-23 21:13:55

客气了，我没做过二次开发，实在帮不上太多。

cainiao07 · 发表于 2014-3-24 22:56:26

chuxf 发表于 2014-3-23 21:13
客气了，我没做过二次开发，实在帮不上太多。

原因已经找到了。可能跟LS-DYAN 主程序有关。主要错误在：
1、if-elseif-endif判断语句未执行，导致循环语句与weibull分布也不能执行
2、LS-DYNA中log函数不可用。原因不明，编译时说是外部函数不可用，我晕。
3、由于有限元计算步之间时间太小，number_seed()函数使用效果不太明显，随机数的随机性还是不太强。我试着增加一个do循环，看能否加大时间步。不过计算机的计算时间暴增了2000倍 TT

楼主你觉得的哥原因里头我是否要申明log函数，或者是不是动态链接库中缺少了这个函数呢，谢谢。

cainiao07 · 发表于 2014-3-25 12:58:56

已解决，在ls-dyna中需要使用log10（）来表示log函数

[数值问题] 程序计算时，参数为什么老是为0呢？

新人勋章

爱心勋章

水王勋章

元老勋章

热心勋章