📄 optexp2.f
字号:
* initial moment22 micron^4/g solvent
moment22(1)=1.801D11
* initial moment30 micron^3/g solvent
moment30(1)=9.203D8
* initial moment31 micron^4/g solvent
moment31(1)=1.801D11
* initial relative supersaturation
relsatn(1)=(concentration(1)-Csat(Temp(0.0D0)))/
& Csat(Temp(0.0D0))
* Initial conditions for derivatives wrt to theta
DO 163 I = 1 , NU
DO 164 J = 1 , NTHETA
derivtheta(1,I,J)=0.0D0
164 CONTINUE
163 CONTINUE
*Noise estimates
conc_variance=3.53D-4
mu00_variance=(0.2D0*moment00(1))**2.0D0
mu10_variance=(0.2D0*seed_moment10(1))**2.0D0
mu01_variance=(0.2D0*seed_moment01(1))**2.0D0
mu11_variance=(0.2D0*seed_moment11(1))**2.0D0
mu20_variance=(0.2D0*seed_moment20(1))**2.0D0
mu02_variance=(0.2D0*seed_moment02(1))**2.0D0
mu21_variance=(0.2D0*seed_moment21(1))**2.0D0
mu12_variance=(0.2D0*seed_moment12(1))**2.0D0
mu22_variance=(0.2D0*moment22(1))**2.0D0
mu30_variance=(0.2D0*moment30(1))**2.0D0
mu31_variance=(0.2D0*moment31(1))**2.0D0
*Simulation parameters
*********************************************************
*
mf=222
itask=1
istate =1
iopt=1
lrw=5000
liw=900
rwork(5) = 1.1D-14
rwork(6) = 1.0D-1
rwork(7) = 1.0D-14
iwork(6) = 100000
rtol=1.0D-6
atol=1.0D-6
itol=1
******************************************************
*
time(1)=0.0D0
T=0.0D0
DO 1200 I=1,(NSTEP-1)
temperature(I)=Temp(T)
Y(1)=moment00(I)
Y(2)=moment10(I)
Y(3)=moment01(I)
Y(4)=moment11(I)
Y(5)=moment20(I)
Y(6)=moment02(I)
Y(7)=moment21(I)
Y(8)=moment12(I)
Y(9)=moment22(I)
Y(10)=moment30(I)
Y(11)=moment31(I)
Y(12)=concentration(I)
Y(13)=seed_moment10(I)
Y(14)=seed_moment01(I)
Y(15)=seed_moment11(I)
Y(16)=seed_moment20(I)
Y(17)=seed_moment02(I)
Y(18)=seed_moment21(I)
Y(19)=seed_moment12(I)
K = 20
DO 263 M = 1 , NU
DO 264 N = 1 , NTHETA
Y(K)=derivtheta(I,M,N)
K = K + 1
264 CONTINUE
263 CONTINUE
CALL lsodes ( MOMENTS,NEQ,Y,T,T+1.0D0,itol,rtol,
& atol,itask,istate,iopt,
& rwork,lrw,iwork,liw, MOMENTSJ, mf )
* print*, 'istate=', istate
* if(istate.EQ.-1) then
* istate=1
* else if (istate.LE.0) then
* print*, "istate=", istate
*
* end if
moment00(I+1)=Y(1)
moment10(I+1)=Y(2)
moment01(I+1)=Y(3)
moment11(I+1)=Y(4)
moment20(I+1)=Y(5)
moment02(I+1)=Y(6)
moment21(I+1)=Y(7)
moment12(I+1)=Y(8)
moment22(I+1)=Y(9)
moment30(I+1)=Y(10)
moment31(I+1)=Y(11)
concentration(I+1)=Y(12)
* concentration_measured(I+1)=concentration(I+1)+
* & betaconc*DRNNOF()
seed_moment10(I+1)=Y(13)
seed_moment01(I+1)=Y(14)
seed_moment11(I+1)=Y(15)
seed_moment20(I+1)=Y(16)
seed_moment02(I+1)=Y(17)
seed_moment21(I+1)=Y(18)
seed_moment12(I+1)=Y(19)
relsatn(I+1)=(concentration(I+1)-Csat(Temp(T)))/
& Csat(Temp(T))
time(I+1)=T
K = 20
DO 363 M = 1 , NU
DO 364 N = 1 , NTHETA
derivtheta(I+1,M,N)=Y(K)
F(Nu*(I-1)+M,N)=Y(k)
K = K + 1
364 CONTINUE
363 CONTINUE
1200 CONTINUE
DO 1370 I = 1, NTHETA
DO 1375 J = 1, NU*(NSTEP-1), Nu
FTV(I,J)=F(J,I)/mu00_variance
FTV(I,J+1)=F(J+1,I)/mu10_variance
FTV(I,J+2)=F(J+2,I)/mu01_variance
FTV(I,J+3)=F(J+3,I)/mu11_variance
FTV(I,J+4)=F(J+4,I)/mu20_variance
FTV(I,J+5)=F(J+5,I)/mu02_variance
FTV(I,J+6)=F(J+6,I)/mu21_variance
FTV(I,J+7)=F(J+7,I)/mu12_variance
FTV(I,J+8)=F(J+8,I)/mu22_variance
FTV(I,J+9)=F(J+9,I)/mu30_variance
FTV(I,J+10)=F(J+10,I)/mu31_variance
FTV(I,J+11)=F(J+11,I)/conc_variance
1375 CONTINUE
1370 CONTINUE
DO 1400 I = 1, NTHETA
DO 1500 J = 1, NTHETA
FTVF(I,J) = 0.0D0
DO 1600 K = 1, Nu*(NSTEP-1)
FTVF(I,J)=FTVF(I,J)+FTV(I,K)*F(K,J)
1600 CONTINUE
1500 CONTINUE
1400 CONTINUE
* DO I=1,NTHETA
* WRITE(*,1650) FTVF(I,1), FTVF(I,2), FTVF(I,3), FTVF(I,4),
* & FTVF(I,5), FTVF(I,6)
* ENDDO
*1650 FORMAT(4(E13.6, 1X))
* PAUSE
*Calculate determinant of FTVF,
*which is symmetric with dimensions NTHETA by NTHETA
DO 200 I = 1, NTHETA
DO 300 J = 1, NTHETA
FAC(I,J) = 0.0D0
300 CONTINUE
200 CONTINUE
CALL DLFTRG (P, FTVF, NTHETA, FAC, LDFAC, IPVT)
CALL DLFDRG (P,FAC,LDFAC,IPVT,DET1,DET2)
detFTVF = DET1*10**(DET2)
fj =-DLOG(detFTVF)
* print*, "fj=", fj
RETURN
END
*****************************************************************
SUBROUTINE cntr(Nvar,jjj,x,gj)
*
* cntr is the subroutine for the constraints. The
* constraints are listed in the following order:
* Nonlinear inequality constraints, linear inequality
* constraints, nonlinear equality constraints, and linear
* equality constraints.
*
* input: Nvar - number of parameters
* jjj - indicates the jjj_th constraint
* x - Nvar-dimensional vector of parameters
*
* output: gj - the jjj_th constraint
*
INTEGER Nvar,jjj, Ntemp2
PARAMETER (Ntemp2 = 8)
REAL*8 x(*),gj, Coeff(Ntemp2+4), umin, umax, Temp
COMMON Coeff
DO 628 I = 1, Nvar
Coeff(I)=x(I)
628 CONTINUE
* Maximum temperature allowed
umax=32.3D0
* Minimum temperature allowed
umin=22.0D0
gj=umin-Temp(160.0D0)
* IF(jjj.LE.Ntemp2) THEN
* gj=Temp(160.0D0/DFLOAT(Ntemp2)*DFLOAT(jjj))-umax
* ELSE
* gj=umin-Temp(160.0D0/DFLOAT(Ntemp2)*DFLOAT(jjj-Ntemp2))
* ENDIF
RETURN
END
****************************************************************
REAL*8 FUNCTION Temp(time)
* Crystallizer temperature setpoint profile for the
* simulation of a batch cooling crystallizer
*
* input: time - minutes
* output: Temp - temperature in degrees Centigrade
INTEGER I, J, Ntemp3
PARAMETER(Ntemp3 = 8)
REAL*8 time, interval, sum
REAL*8 Coeff(Ntemp3+4)
COMMON Coeff
interval = DFLOAT(160/Ntemp3)
sum=0.0D0
DO 2666 I = 1, (Ntemp3-1)
IF(time.LE.(DFLOAT(I)*interval)) THEN
DO 2667 J = 1, (I-1)
sum=sum+Coeff(J)
2667 CONTINUE
Temp=32.0D0+sum*interval+
& Coeff(I)*(time-DFLOAT(I-1)*interval)
RETURN
ENDIF
2666 CONTINUE
DO 2668 I = 1, (Ntemp3-1)
sum=sum+Coeff(I)
2668 CONTINUE
Temp=32.0D0+sum*interval+Coeff(Ntemp3)*
& (time-DFLOAT(Ntemp3-1)*interval)
RETURN
END
****************************************************************
SUBROUTINE MOMENTS(NEQ,T,Y,Yprime)
INTEGER NEQ, I, J, K
INTEGER NTHETA, NU
REAL*8 T, Y(NEQ), Yprime(NEQ)
REAL*8 r0, alpha, mu00, kg1, g1, g2, kg2, kb, b
REAL*8 Csat, birth, birth2
REAL*8 growth1, growth2, growth3, growth4, Temp
REAL*8 JACOBIAN(12,12), W(12,6)
REAL*8 WPRIME(12,6), F_THETA(12,6)
COMMON /EXP_DATA/r0, alpha, mu00
COMMON /GROWTH_DATA1/kg1, g1
COMMON /GROWTH_DATA2/kg2, g2
COMMON /BIRTH_DATA/kb, b
NTHETA=6
NU=12
*Moments
* 00 moment, number of particles/g solvent/minute
Yprime(1)=birth(Y(12),Csat(Temp(T)),Y(7))
* 10 moment, microns/g solvent/minute
Yprime(2)=growth1(Y(12),Csat(Temp(T)))*Y(1)
* 01 moment, microns/g solvent/minute
Yprime(3)=growth2(Y(12),Csat(Temp(T)))*Y(1)
* 11 moment, microns^2/g solvent/minute
Yprime(4)=growth1(Y(12),Csat(Temp(T)))*Y(3)
& + growth2(Y(12),Csat(Temp(T)))*Y(2)
* 20 moment, microns^2/g solvent/minute
Yprime(5)=2*growth1(Y(12),Csat(Temp(T)))*Y(2)
* 02 moment, microns^2/g solvent/minute
Yprime(6)=2*growth2(Y(12),Csat(Temp(T)))*Y(3)
* 21 moment, microns^3/g solvent/minute
Yprime(7)=2*growth1(Y(12),Csat(Temp(T)))*Y(4)+
& growth2(Y(12),Csat(Temp(T)))*Y(5)
* 12 moment, microns^3/g solvent/minute
Yprime(8)=growth1(Y(12),Csat(Temp(T)))*Y(6)+
& 2*growth2(Y(12),Csat(Temp(T)))*Y(4)
* 22 moment, microns^4/g solvent/minute
Yprime(9)=2*growth1(Y(12),Csat(Temp(T)))*Y(8)+
& 2*growth2(Y(12),Csat(Temp(T)))*Y(7)
* 30 moment, microns^4/g solvent/minute
Yprime(10)=3*growth1(Y(12),Csat(Temp(T)))*Y(5)
* 31 moment, microns^4/g solvent/minute
Yprime(11)=3*growth1(Y(12),Csat(Temp(T)))*Y(7)+
& growth2(Y(12),Csat(Temp(T)))*Y(10)
*Mass balance
* d(conc(t))/dt, g solute/g solvent/minute
Yprime(12)=-alpha*2.0D0*growth1( Y(12), Csat(Temp(T)) )
& *Y(4) - alpha*growth2( Y(12),Csat(Temp(T) ) )*Y(5)
*Moments for seed crystals only
* 10 moment, microns/g solvent/minute
Yprime(13)=growth1(Y(12),Csat(Temp(T)))*mu00
* 01 moment, microns/g solvent/minute
Yprime(14)=growth2(Y(12),Csat(Temp(T)))*mu00
* 11 moment, microns^2/g solvent/minute
Yprime(15)=growth1(Y(12),Csat(Temp(T)))*Y(14)
& + growth2(Y(12),Csat(Temp(T)))*Y(13)
* 20 moment, microns^2/g solvent/minute
Yprime(16)=2*growth1(Y(12),Csat(Temp(T)))*Y(13)
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -