Я полный новичок в Фортране и работаю над кодом, который решает кинетический механизм, решая дифференциальные уравнения на разных временных шагах с использованием решателя дифференциальных уравнений.Вот ссылка для загрузки zip-файла со всем проектом:
http://www.filedropper.com/fortranmicropyrolyzersetup
Входные переменные для дифференциального решателя определены следующим образом в коде:
! Declaration of variables
implicit none
EXTERNAL :: FEXSB_AUTO, JEX_SB
integer :: neq,Mf,lrw,liw,iwork,itol,itask,istate,iopt !solver parameters
integer :: j,jjk,m !Counters
double precision :: ATOL,RTOL, RWORK !solver parameters
double precision :: T, TOUT !starting time (s), timestep time (s)
double precision :: Y, w !molar fraction of biomass (-), mass fraction of biomass (-)
double precision :: y_gas, w_gas, Conc !molar fraction in gas phase (-), concentration in gas phase (mol/m³), mass fraction in gas phase (-)
double precision :: n(speciescount) !number of moles (used temporarily to calculate initial molar fracions)
character*5 :: simulationNumber
!Setting the solver parameters
neq = SpeciesCount !number of equations
ITOL = 1 !RTOL and ATOL are integers
RTOL = 1.0D-8 !Relative tolerance
ATOL = 1.0D-15 !Absolute tolerance
ITASK = 1
ISTATE = 1
LRW = 22 + 9*SpeciesCount + 2*SpeciesCount**2 !Array sizing (see VODE)
LIW = 30+SpeciesCount !Array sizint (see VODE)
MF = 22 !Use BDF with finite difference Jacobian
IOPT = 1 !Optional input specified
Iwork(6) = 7000 !Increase maximum iteration steps from 500 to 2000 otherwise the solver does not converge
Затем дифференциальный решатель вызывается в цикле do для каждого временного шага следующим образом:
! Solve reactor equations (see FEXSB) and advance time, until stop criterium is met
TimeStep = 0 ! dimensionless time step
!DO while(wm(1).lt.0.9999) -> previously used stop criterium
DO while((y_gas(1).lt.0.99999).OR.(TimeStep.lt.10)) !stop criterium: molar fraction Helium = 0.9999 AND do at least 100 timesteps
TimeStep = TimeStep + 1
write(*,*)"Doing iteration for time step",TimeStep
call DVODE(FEXSB_AUTO,NEQ,Y,T,TOUT,ITOL,RTOL,ATOL,ITASK,ISTATE,IOPT,RWORK,LRW,IWORK,LIW,JEX_SB,MF) !solve reactor equations to Y
!CALL DVODE(FEXSB_AUTO,NEQ,Y,T,TOUT,ITOL,RTOL,ATOL,ITASK,ISTATE,IOPT,RWORK,LRW,IWORK,LIW,JEX_SB,MF,RPAR,IPAR)
! calculate w, y_gas, w_gas and Conc from Y
do j = 1, SpeciesCount
if (Y(j) .lt. 1.0D-10) then ! round to zero below 10e-10 to avoid negative numbers and numerical problems with small numbers
Y(j) = 0.0D0
endif
if (MolarFlowRate(j) .lt. 1.0D-10) then
MolarFlowRate(j) = 0.0D0
endif
w(j) = Y(j)*n0_tot*amms(j) / (mass_sample)
y_gas(j) = MolarFlowRate(j) / TotalMolarFlowRate
w_gas(j) = MolarFlowRate(j)*amms(j) / (1000*MassFlowRate) ! factor 1000 to put molar mass in kg/mol instead of g/mol
Conc(j) = MolarFlowRate(j) / VolumetricFlowRate
end do
Средство вычисления дифференциальных уравнений не успешно завершает цикл do строки:
DO while((y_gas(1).lt.0.99999).OR.(TimeStep.lt.10))
Кажется, что переменные распознаются для первого временного шага при вызове ODE Solver в строке:
call DVODE(FEXSB_AUTO,NEQ,Y,T,TOUT,ITOL,RTOL,ATOL,ITASK,ISTATE,IOPT,RWORK,LRW,IWORK,LIW,JEX_SB,MF)
И решатель успешно завершает работупервая итерация.После увеличения временного шага функция вызова снова работает, но на этот раз переменные как-то не распознаются.Когда я остановил код для отладки того, что не так после первого временного шага, я понял, что переменные, необходимые для DVODE, больше не имеют заданных значений, так или иначе они удаляются после первой успешной итерации.
Что может бытьвызывая эту проблему?
*DECK DVODE
SUBROUTINE DVODE (F, NEQ, Y, T, TOUT, ITOL, RTOL, ATOL, ITASK,
1 ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JAC, MF,
2 RPAR, IPAR)
EXTERNAL F, JAC
DOUBLE PRECISION Y, T, TOUT, RTOL, ATOL, RWORK, RPAR
INTEGER NEQ, ITOL, ITASK, ISTATE, IOPT, LRW, IWORK, LIW,
1 MF, IPAR
DIMENSION Y(*), RTOL(*), ATOL(*), RWORK(LRW), IWORK(LIW),
1 RPAR(*), IPAR(*)
!-----------------------------------------------------------------------
c dvode: Variable-coefficient Ordinary Differential Equation solver,
! with fixed-leading-coefficient implementation.
! This version is in double precision.
!
! DVODE solves the initial value problem for stiff or nonstiff
! systems of first order ODEs,
! dy/dt = f(t,y) , or, in component form,
! dy(i)/dt = f(i) = f(i,t,y(1),y(2),...,y(NEQ)) (i = 1,...,NEQ).
! DVODE is a package based on the EPISODE and EPISODEB packages, and
! on the ODEPACK user interface standard, with minor modifications.
!-----------------------------------------------------------------------
! Authors:
! Peter N. Brown and Alan !. Hindmarsh
! Center for Applied Scientific Computing, L-561
! Lawrence Livermore National Laboratory
! Livermore, CA 94551
! and
! George D. Byrne
! Illinois Institute of Technology
! Chicago, IL 60616
!-----------------------------------------------------------------------
Любая помощь будет принята с благодарностью.Обратите внимание, что я новичок в Фортране.Если я предоставлю какую-либо дополнительную информацию, чтобы помочь вам ответить на мой вопрос, пожалуйста, не стесняйтесь, дайте мне знать.Заранее спасибо!