在 ifort 中标记以避免将不一致的数据类型传递给 Fortran 子例程

问题描述

我正在 Fortran 90 中编写第 3 方代码，此代码过去运行良好。然后我工作的集群最近有更新，我不确定细节（我认为编译器保持不变），然后这个 3rd 方代码停止工作。阅读代码我看到一些整数值，如 1 和 -1 被传递给没有变量的子例程。示例：

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt 
import math as math
import pandas as pd

#data set 1
t_data1 = [0,0.08333,0.5,1,4,22.61667]
x_data1 = [0,4.91e-5,4.57e-5,4.74e-5,4.17e-5,2.76e-5]
x_data1mgl = [0,3.48,3.24,3.36,2.96,1.96]

#data set 2
t_data2 = [0,22.8167]
x_data2 = [0,5.92e-5,5.7e-5,5.64e-5,5.30e-5,4.6e-5]
x_data2mgl = [0,4.2,4.04,3.76,2.88]
      
#combine data and time in the same dataframe
data1 = pd.DataFrame({'time':t_data1,'x1':x_data1})
data2 = pd.DataFrame({'time':t_data2,'x2':x_data2})
data1.set_index('time',inplace=True)
data2.set_index('time',inplace=True)


# merge dataframes
data = data1.join(data4,how='outer')

# simulation time points
dftp = pd.DataFrame({'time':np.linspace(0,25,51)})
dftp.set_index('time',inplace=True)

# merge dataframes
df = data.join(dftp,how='outer')

# get True (1) or False (0) for measurement
#df['meas'] = (df['x'].values==df['x'].values).astype(int)
z1 = (df['x1']==df['x1']).astype(int)
z2 = (df['x2']==df['x2']).astype(int)


# replace NaN with zeros
df0 = df.fillna(value=0)

#Estimator Model
m = GEKKO()#remote=True)
#m = GEKKO(remote=False) # remote=False to produce local folder with results

m.time = df.index.values

# measurements - Gekko arrays to store measured data
xm = m.Array(m.Param,2)
xm[0].value = df0['x1'].values
xm[1].value = df0['x2'].values

hocl_init_val=m.Array(m.Param,2)
hocl_init_val[0].value= 8.01e-5
hocl_init_val[1].value= 8.01e-5

nh3_init_val=m.Array(m.Param,2)
nh3_init_val[0].value= 1.37e-3 
nh3_init_val[1].value= 6.82e-4

h_init_val=m.Array(m.Param,2)
h_init_val[0].value= 6.31e-8 
h_init_val[1].value= 2e-8

#h2co3_init_val=m.Array(m.Param,2)
#h2co3_init_val[0].value= 5.39e-4
#h2co3_init_val[1].value= 1.88e-4

hco3_init_val=m.Array(m.Param,2)
hco3_init_val[0].value= 3.46e-3
hco3_init_val[1].value= 3.80e-3

#co32_init_val=m.Array(m.Param,2)
#co32_init_val[0].value= 2.16e-6
#co32_init_val[1].value= 7.5e-6

alk_init_val=m.Array(m.Param,2)
alk_init_val[0].value= 3.46e-3
alk_init_val[1].value= 3.82e-3

#adjustable parameters

x1=3.1684e-5
x2=3.1308e-5

DOC10=m.Array(m.FV,2)
DOC10[0].value= x1
DOC10[1].value= x1*0.5

DOC20=m.Array(m.FV,2)
DOC20[0].value= x2
DOC20[1].value= x2*0.5


kdoc1 = m.FV(1e6,lb=0.01,ub=1e7) #m.Const(1.2334e6)
kdoc2 = m.FV(1e9,ub=1e10) #m.Const(3.8809e9)
#DOC10 = m.FV(1e-5,lb=1e-7,ub=1e-2)
#DOC20 = m.FV(1e-5,ub=1e-2)


#variables initialization
# hocl=m.Array(m.Var,2,lb=0,ub=1e-4)
# nh3=m.Array(m.Var,ub=1.5e-4)
# nh2cl=m.Array(m.Var,ub=1e-4)
# C2m=m.Array(m.Var,2)
# cC2=m.Array(m.Var,2)
# nhcl2=m.Array(m.Var,2)
# h=m.Array(m.Param,2)
# #oh=m.Array(m.Var,2)
# I=m.Array(m.Var,2)
# ocl=m.Array(m.Var,2)
# nh4=m.Array(m.Var,2)
# #h2co3=m.Array(m.Var,2)
# hco3=m.Array(m.Var,2)
# #co32=m.Array(m.Var,2)
# alk=m.Array(m.Var,2)
# DOC1=m.Array(m.SV,2)
# DOC2=m.Array(m.SV,2)

#to support intermediate
co32=[None]*2
h2co3=[None]*2

oh=[None]*2

r1=[None]*2
r2=[None]*2
r3=[None]*2
r4=[None]*2
r5=[None]*2
r6=[None]*2
r7=[None]*2
r8=[None]*2
r9=[None]*2
r10=[None]*2
r11=[None]*2
r12=[None]*2
r13=[None]*2
r14=[None]*2
r15=[None]*2
r16=[None]*2

#Define GEKKO variables that determine if time 
#point contains data to be used in regression
zm = m.Array(m.Param,2)
zm[0].value=z1
zm[1].value=z2

# fit to measurement
x=m.Array(m.Param,2) 
x[0].value=x_data1[0]
x[1].value=x_data4[0]

k5=m.Array(m.Var,2)

k1 = m.Const(1.5e10)
k2 = m.Const(7.6e-2)
k3 = m.Const(1e6)
k4 = m.Const(2.3e-3)
k6 = m.Const(2.2e8)
k7 = m.Const(4e5)
k8 = m.Const(1e8)
k9 = m.Const(3e7)
k10 = m.Const(55)
k11 = m.Const(3.16e-8*1e10)
k12 = m.Const(1e10)
k13 = m.Const(5.01e-10*1e10)
k14 = m.Const(1e10)

for i in range(2):

    hocl[i] = m.Var(value=hocl_init_val[i],ub=1e-2)
    nh3[i] = m.Var(value=nh3_init_val[i],ub=1.5e-3)
    nh2cl[i] = m.Var(value=x[i],ub=3e-4)
    
    C2m[i] = m.Param(xm[i],lb=0)
    meas = m.Param(zm[i])
    m.Minimize(meas*((nh2cl[i]-C2m[i]))**2)
    
    nhcl2[i] = m.Var(value=0)
    
    h[i] = m.Param(value=h_init_val[i])
    
    I[i] = m.Var(value=0)
    ocl[i]= m.Var(value=0)
    nh4[i] = m.Var(value=0)
    #h2co3[i] = m.Var(value=h2co3_init_val[i])
    hco3[i] = m.Var(value=hco3_init_val[i])
    #co32[i] = m.Var(value=co32_init_val[i])
    alk[i] = m.Var(value=alk_init_val[i])
    DOC1[i] = m.SV(value=DOC10[i])#,fixed_initial=False)
    DOC2[i] = m.SV(value=DOC20[i])#,fixed_initial=False)
    cC2[i] = m.Var(value=0)
    
    oh[i] = m.Intermediate(1e-14/h_init_val[i])

    r1[i] = m.Intermediate(k1 * hocl[i] * nh3[i])
    r2[i] = m.Intermediate(k2 * nh2cl[i])
    r3[i] = m.Intermediate(k3 * hocl[i] * nh2cl[i])
    r4[i] = m.Intermediate(k4 * nh3[i])
    r5[i] = m.Intermediate(k5[i] * nh2cl[i] * nh2cl[i])
    r6[i] = m.Intermediate(k6 * nhcl2[i] * nh3[i]* h[i])
    r7[i] = m.Intermediate(k7 * nhcl2[i] * oh[i])
    r8[i] = m.Intermediate(k8 * I[i] * nhcl2[i])
    r9[i] = m.Intermediate(k9 * I[i] * nh2cl[i])
    r10[i] = m.Intermediate(k10 * nh2cl[i] * nhcl2[i])
    
    r11[i] = m.Intermediate(k11*hocl[i])
    r12[i] = m.Intermediate(k12*h[i]*ocl[i])
    r13[i] = m.Intermediate(k13*nh4[i])
    r14[i] = m.Intermediate(k14*h[i]*nh3[i])
    
    r15[i] = m.Intermediate(kdoc1*DOC1[i]*nh2cl[i])
    r16[i] = m.Intermediate(kdoc2*DOC2[i]*hocl[i])

    co32[i] = m.Intermediate(5.01e-11*hco3[i]/h[i])
    h2co3[i] = m.Intermediate(hco3[i]*h[i]/5.01e-7)
    
    t = m.Param(value=m.time)
    m.Equation(hocl[i].dt()== -r1[i] + r2[i] - r3[i] + r4[i] + r8[i] - r11[i] + r12[i] - r16[i])
    m.Equation(nh3[i].dt()== -r1[i] + r2[i] + r5[i] - r6[i] + r13[i] - r14[i])# + r11[i])
    m.Equation(nh2cl[i].dt()== r1[i] - r2[i] - r3[i] + r4[i] - r5[i] + r6[i] - r9[i] - r10[i] - r15[i])
    m.Equation(nhcl2[i].dt()== r3[i] - r4[i] + r5[i] - r6[i] - r7[i] - r8[i] - r10[i])
    #m.Equation(h[i].dt()== 0)
    #m.Equation(oh[i] == 1e-14/h[i])
    m.Equation(I[i].dt()== r7[i] - r8[i] - r9[i])
    #m.Equation(ocl[i] == (3.16e-8*hocl[i])/h[i])
    #m.Equation(nh4[i] == nh3[i]*h[i]/5.01e-10)
    m.Equation(ocl[i].dt()==r11[i]-r12[i])
    m.Equation(nh4[i].dt()==-r13[i]+r14[i])
    #m.Equation(co32[i] == (5.01e-11*hco3[i])/h[i])
    #m.Equation(h2co3[i] == (hco3[i])*h[i]/5.01e-7)
    m.Equation(hco3[i] == alk[i] - 2*5.01e-11*hco3[i]/h[i] - oh[i] + h[i])
    m.Equation(alk[i].dt()== 0)
    
    m.Equation(DOC1[i].dt()== -r15[i])
    m.Equation(DOC2[i].dt()== -r16[i])
    
    m.Equation(cC2[i] == 51500*nh2cl[i])
    
    m.Equation(k5[i] == 2.5e7*h[i]+4e4*(hco3[i]*h[i]/5.01e-7)+800*hco3[i])

    #m.Connection(DOC1[i],DOC10[i],pos1=1,pos2=1,node1=1,node2=1)
    #m.Connection(DOC2[i],DOC20[i],node2=1)
    
#Application options
m.options.soLVER = 3 #IPOPT solver
m.options.EV_TYPE = 2 #absolute error
#m.options.NODES = 3 #collocation nodes (2,5)
#m.options.RTOL = 1E-6

if True:
    kdoc1.STATUS=1
    kdoc2.STATUS=1
    DOC10[i].STATUS=1
    DOC20[i].STATUS=1

#m.options.IMODE = 7
#m.options.RTOL = 1E-6
#m.solve(dis=True)

m.options.IMODE = 5 #Dynamic Simultaneous - estimation = MHE
m.options.COLDSTART=2
m.options.TIME_SHIFT=0

m.open_folder()
m.solve(disp=True)

print('Final SSE Objective: ' + str(m.options.objfcnval))

print('Solution')
print('kdoc1 = ' + str(kdoc1.value[0]))
print('kdoc2 = ' + str(kdoc2.value[0]))
print('DOC10_w1 = ' + str(DOC10[0].value[0]))
print('DOC20_w1 = ' + str(DOC20[0].value[0]))
print('DOC10_w2 = ' + str(DOC10[1].value[1]))
print('DOC20_w2 = ' + str(DOC20[1].value[1]))

plt.figure(1,figsize=(8,5))
plt.subplot(2,1)
plt.plot(m.time,nh2cl[0],'b.--',label='Predicted 1')
plt.plot(m.time,nh2cl[1],'r.--',label='Predicted 2')

plt.plot(t_data1,x_data1,'bx',label='Meas 1')
plt.plot(t_data2,x_data2,'rx',label='Meas 2')
#plt.legend(loc='best')
plt.ylabel('mol/L')
plt.legend(loc='center right',bBox_to_anchor=(1.45,0.5),ncol=2) #shadow=True,plt.subplot(2,2)
plt.plot(m.time,cC2[0].value,'b',label ='C2_M1')
plt.plot(m.time,cC2[1].value,'r',label ='C2_M2')
plt.plot(t_data1,x_data1mgl,x_data2mgl,label='Meas 2')

plt.ylabel('mg Cl2/L')
plt.xlabel('time (h)')
plt.legend(loc='center right',bBox_to_anchor=(1.4,plt.figure(2,figsize=(12,8))
plt.subplot(4,3,hocl[i].value,label ='hocl')
plt.legend()

plt.subplot(4,nh3[i].value,label ='nh3')
plt.legend()

plt.subplot(4,3)
plt.plot(m.time,nhcl2[i].value,label ='nhcl2')
plt.legend()

plt.subplot(4,4)
plt.plot(m.time,h[i].value,label ='h')
plt.legend()

plt.subplot(4,5)
plt.plot(m.time,oh[i].value,label ='oh')
plt.legend()

plt.subplot(4,6)
plt.plot(m.time,I[i].value,label ='I')
plt.legend()

plt.subplot(4,7)
plt.plot(m.time,ocl[i].value,label ='ocl')
plt.legend()

plt.subplot(4,8)
plt.plot(m.time,nh4[i].value,label ='nh4')
plt.legend()

plt.subplot(4,9)
plt.plot(m.time,h2co3[i].value,label ='h2co3')
plt.legend()
plt.xlabel('time (h)')

plt.subplot(4,10)
plt.plot(m.time,hco3[i].value,label ='hco3')
plt.legend()
plt.xlabel('time (h)')

plt.subplot(4,11)
plt.plot(m.time,co32[i].value,label ='co32')
plt.legend()
plt.xlabel('time (h)')

plt.subplot(4,12)
plt.plot(m.time,alk[i].value,label ='alk')
plt.legend()
plt.xlabel('time (h)')

plt.show()

以及稍后在主文件中：

integer,parameter :: sp = selected_int_kind(18)

subroutine sub(var)

integer(kind=sp),intent(in):: var

<do something>

return
end subroutine

用 call sub(1) 检查变量我已经看到在 gdb sub 里面不是 1，而是一些大整数。这是由于整数类型不一致造成的。

我知道如何解决这个问题，var 会发生这种情况。

我有两个问题：

• 为什么现在集群更新后会出现这个问题？

• ifort 编译器中是否有一个标志可以在传递不一致的类型时返回错误或警告？

解决方法

有一个 ifort 标志可以提供帮助 - -warn interface。如果启用此选项，编译器将为不在模块中的外部过程创建一个接口块（经典 F77 样式）。然后，当对过程进行调用时，如果没有显式接口，编译器将查看是否有生成的接口并比较它们，如果它们不匹配，则向您发出警告。它并不完美 - 它取决于首先编译的被调用例程。

至于为什么会发生，这完全取决于传递参数后内存中发生了什么。

最好的解决方案是为所有过程提供一个明确的接口。当前版本支持 Fortran 2018 的 IMPLICIT NONE(EXTERNAL)，它强制您为您调用的所有过程具有显式接口（或 EXTERNAL 声明）。还有一个编译选项可以打开它，-warn externals，但在版本 18 中都没有。