cgxbfgr.for
c
cFile-name PHOTWO.htm
c
SUBROUTINE PHOTWO (LOUT,KK,NATJ,NMAX,JJ,LLNRGY,LENTWP,LEVELA,
1 ACTIVE,ICKWRK,RCKWRK,VOL,WT,
2 ATIM,RTIM,ATOL,RTOL,ABSOL,RELAT,
3 ABOVE,BELOW,BUFFER,TWPWK,IPVT,SCRTCH,
4 S,SN,RES,RESN,A,SU,AP,H0K,UFAC,DFAC,
6 P,T,NINIT,NJAC,ITJAC,IRETIR,NUMDT,DTTP,
7 DTMIN,DTMAX,QDOT,ICALL1,T0,DBGA,LERR,TPSEM)
C
C*****precision > double
IMPLICIT DOUBLE PRECISION (A-H, O-Z), INTEGER (I-N)
REAL APP,SUU
C*****END precision > double
C*****precision > single
C IMPLICIT REAL (A-H,O-Z), INTEGER (I-N)
C*****END precision > single
C
DIMENSION LEVEL(2), ICKWRK(*), RCKWRK(*), WT(*), ABOVE(*),
1 BELOW(*), BUFFER(*), TWPWK(*), IPVT(*), SCRTCH(KK,*),
2 ITWPWK(3), S(*), SN(*), RES(*), RESN(*), A(NATJ,*),
3 SU(*), AP(*), H0K(*), LEVELA(2)
DIMENSION X(NMAX)
C
INTEGER CALL, CALLS, CALL1
C
LOGICAL ACTIVE(*), MARK(1), LENRGY,
1 ERROR, FUNCTN, JACOBN, REENTR, SOLVE, STORE, SUCCES,
2 ADAPT, SHOW, SAVE, UPDATE, ENERGY, LTIME, DBG, DBGA,
3 FRSTRS, LLNRGY, LERR, TPSEM
C
CHARACTER*6 NAMSUB,NAMFUN
COMMON /NAMFN/NAMSUB,NAMFUN
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C INPUT-
C LOUT - unit for printed output
C KK - number of chemical species.
C NATJ - number of dependent variables solved. NATJ=KK if
C temperature is solved, KK-1 othewise.
C JJ - dummy (no. of cells for future use).
C LLNRGY - if LLNRGY=.TRUE. then the energy equation is to be
C solved, otherwise a fixed temperature is used.
C LENTWP - size of real workspace for TWOPNT.
C LEVEL - LEVEL(1) controls the printed output from TWOPNT and
C LEVEL(2) controls the printing from photwo. LEVEL(1) must
C always be greater than or equal to LEVEL(2). LEVEL may
C have values of 0, 1, 2, or 3.
C dimension LEVEL(2).
C ACTIVE - active(*) is a logical array that controls which variables
C TWOPNT will attempt adaptive meshing. currently,
C there is no mesh, therefore all active(*)=.false. (in
C future a mesh will be possible) logical active(kk).
C VOL - volume of the current cell (expand to an array later)
C WT - the array of species molecular weights.
C cgs units - gm/mole
C dimension wt(*) at least kk.
C ATIM - absolute convergence criteria for the newton iteration
C as used for the time steps.
C RTIM - relative convergence criteria for the newton iteration
C as used for the time steps.
C ATOL - absolute convergence criteria for the newton iteration.
C RTOL - relative convergence criteria for the newton iteration.
C ABSOL - absolute perturbation for computing jacobian.
C RELAT - relative perturbation for computing jacobian.
C ABOVE - array of upper bounds for the dependent variables. used
C by TWOPNT to control the damping. above has the same
C structure as the solution vector s(*). above(*) should be
C set to values that are above acceptable solution values.
C dimension above(*) at least NATJ.
C BELOW - array of lower bounds for the dependent variables. used
C by TWOPNT to control the damping. below has the same
C structure as the solution vector s(*). below(*) should be
C set to values that are below acceptable solution values.
C dimension below(*) at least NATJ.
C S - dependent variable array. on input the solution estimate
C is in s. the temperature is stored in t=s(nt), and the
C mass fractions are in y(k)=s(nys+k).
C dimension s(*) at least NATJ.
C SU - phoenics calculated residuals (without chemical sources).
C AP - phoenics calculated diagonal terms of the jacobian.
C H0K - reference enthalpies (for calculating the heat source).
C UFAC - factor by which the internal false-timestep may be increased
C each time an increase is appropriate.
C DFAC - factor by which the internal false-timestep may be decreased
C each time an decrease is appropriate.
C P - the pressure.
C cgs units - dynes/cm**2
C T - the temperature. for fixed-temperature problems this
C is the user-specified temperature.
C cgs units - k
C NINIT - the number of time-steps to be taken before the steady
C problem (no internal false time-step) is attempted.
C NJAC - no. of iterations allowed before the jacobian is
C recalculated when solving the steady problem.
C ITJAC - no. of iterations allowed before the jacobian is
C recalculated when solving the transient problem.
C IRETIR - number of steps taken before the time-step may be
C increased.
C NUMDT - the number of time steps to be taken upon failure of a
C newton iteration. numdt is in effect when either the
C temperature-fixed problem or the temperature-varying
C problem is being solved.
C DTTP - the value of the timestep that is taken when the newton
C iteration fails.
C cgs units - sec
C DTMIN - lower limit on the false-timesteps
C DTMAX - upper limit on the false-timesteps
C ICALL1 - first solution stage to attempted; if 1 then solve the
C fixed temperature problem before attempting solution
C of the energy equation, if 2 then go directly to solve
C the energy equation.
C
C WORK AND SCRATCH SPACE-
C RCKWRK - floating point chemkin work space.
C dimensioning - see chemkin documentation.
C ICKWRK - integer chemkin work space.
C dimensioning - see chemkin documentation.
C BUFFER - work space into which TWOPNT writes data and through which
C data is returned to TWOPNT through the reverse
C communication interface. see TWOPNT documentation.
C dimension below(*) at least NATJ.
C TWPWK - work space for TWOPNT.
C dimension wnewtn(lentwp).
C IPVT - array of pivots used by linpack lu factorization and solve
C routines.
C dimension ipvt(*) at least NATJ.
C SCRTCH - scratch space used by the function, jacobian, and print.
C dimension scrtch(kk,5)
C SN - dependent variable array at previous timestep. the
C structure is the same as s.
C dimension sn(*) at least NATJ.
C RES - residuals of the governing equations as evaluated at
C s(n). the residual of the kk species equations is in
C res(nys+k), the energy equation in res(nt,j). if
C inergy=.false. then the energy equation is replaced
C by the given temperature.
C dimension res(*) at least NATJ.
C A - storage space for the jacobian matrix, and its lu factors.
C dimension a(NATJ,NATJ)
C
C OUTPUT-
C S - dependent variable array. on output s(*) contains the
C solution. the temperature is stored in t=s(nt), and the
C mass fractions are in y(k)=s(nys+k)
C dimension s(*) at least NATJ.
C QDOT - the heat release-rate.
C cgs units - ergs/sec
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
DATA GRAD, CURV, ADAFLR /3*1.0/, ADAPT/.FALSE./
C
NAMSUB='PHOTWO'
if(dbga) then
level(1)=levela(1)
level(2)=levela(2)
else
level(1)=0
level(2)=0
endif
dbg=dbga.or.level(2).gt.0
if(dbg) then
call banner(1,'photwo',140422)
app=ap(1); suu=su(1) ! app, suu are local REAL variables
call writ2r('ap ',app,'su ',suu)
level(1)=2
level(2)=2
endif
C
C SET NDT AND VDT FOR THE TEMPERATURE FIXED PROBLEM
C
NDT = NUMDT
VDT = DTTP
LENRGY = LLNRGY
LTIME = .FALSE.
LERR=.FALSE.
C
CALL PHOFUN (KK, NATJ, LENRGY, P, T, VOL, ICKWRK, RCKWRK,
1 WT, SCRTCH(1,1), SCRTCH(1,2), SN, S,
2 SCRTCH(1,3), SU, AP, H0K, QDOT, LTIME, VDT, T0)
C
C DECIDE HOW MANY TIMES TO CALL TWOPNT
C
ENERGY = LENRGY
IF(ENERGY) THEN
CALLS = 2
ELSE
CALLS = 1
ENDIF
CALL1=MIN(CALLS,ICALL1)
C
C TOP OF THE LOOP OVER CALLS TO TWOPNT.
C
DO 1000 CALL = CALL1, CALLS
C
IF(CALL.EQ.2 .AND. DBG)
1 CALL WRIT40(' PHOTWO - Isothermal done; solve energy ')
C
IF(ENERGY) LENRGY = (CALL .NE. 1)
C
REENTR = .FALSE.
C
IPASSS = 1
500 CONTINUE
C
CALL TWOPNT (ERROR, LOUT, LEVEL, 3, ITWPWK, LENTWP, TWPWK,
1 ABOVE, ACTIVE, ADAPT, BELOW, 0, BUFFER, NATJ,
2 CONDIT, FUNCTN, JACOBN, MARK, X, IPASS, IPASSS,
3 NADP, NMAX, JJ, REENTR, IREPRT, SAVE, 0,
4 SHOW, SOLVE, ATOL, NJAC, RTOL, NINIT, NUMDT,
5 IRETIR, STORE, SUCCES, ATIM, ITJAC, DFAC, RTIM,
6 LTIME, UFAC, DTMAX, DTMIN, ADAFLR, GRAD, CURV,
7 VDT, DT, UPDATE, S)
C
C--- Reset initial time-step to min of last time-step and externally
C set value
DTTP = MIN(DTTP, DT)
C
IF(ERROR) THEN
if(dbg) then
write(lout,'('' photwo - error'')')
write(lout,'('' Solution: '')')
write(lout,'(1x,1p6e10.3)') (s(i),i=1,natj)
endif
IF(TPSEM)
1 CALL WRIT40(' PHOTWO - TWOPNT failed to converge ')
LERR=.TRUE.
RETURN
ELSEIF(.NOT.SUCCES .AND. .NOT.REENTR) THEN
if(dbg) then
write(lout,'('' photwo - solution failed'')')
write(lout,'('' Solution: '')')
write(lout,'(1x,1p6e10.3)') (s(i),i=1,natj)
endif
CALL WRIT40(' PHOTWO - TWOPNT failed to converge ')
LERR=.TRUE.
RETURN
ELSEIF(REENTR) THEN
C
IF(FUNCTN) THEN
C
if(dbg) then
write(lout,'('' photwo - residual'')')
endif
if(level(2).gt.1) then
write(lout,'('' Solution (calc. res.): '')')
write(lout,'(1x,1p6e10.3)') (buffer(i),i=1,natj)
write(lout,'('' Ap: '')')
write(lout,'(1x,1p6e10.3)') (ap(i),i=1,natj)
write(lout,'('' Su: '')')
write(lout,'(1x,1p6e10.3)') (su(i),i=1,natj)
endif
CALL PHOFUN (KK, NATJ, LENRGY, P, T, VOL, ICKWRK, RCKWRK,
1 WT, SCRTCH(1,1), SCRTCH(1,2), SN, BUFFER,
2 RES, SU, AP, H0K, QDOT, LTIME, DT, T0)
C
C*****precision > double
CALL DCOPY (NATJ, RES, 1, BUFFER, 1)
cc IF(FRSTRS) CALL DCOPY (NATJ, RES, 1, SCRTCH(1,3), 1)
C*****END precision > double
C*****precision > single
C CALL SCOPY (NATJ, RES, 1, BUFFER, 1)
C IF(FRSTRS) CALL SCOPY (NATJ, RES, 1, SCRTCH(1,3), 1)
C*****END precision > single
FRSTRS=.FALSE.
ELSEIF(JACOBN) THEN
if(dbg) write(lout,'('' photwo - jacobian'')')
CALL PHOJAC (KK, NATJ, LENRGY, ABSOL, RELAT, P, T, VOL,
1 ICKWRK, RCKWRK, WT, SCRTCH, SN, S,
3 RES, RESN, A, SU, AP, H0K, LTIME, DT, T0)
C*****precision > double
CALL DGECO (A, NATJ, NATJ, IPVT, RCOND, RESN)
C*****END precision > double
C*****precision > single
C CALL SGECO (A, NATJ, NATJ, IPVT, RCOND, RESN)
C*****END precision > single
IF(RCOND .LE. 0.0E0) THEN
CALL WRIT40(
1 ' PHOTWO - Singular Jacobian (reduce DTFs)')
LERR=.TRUE.
RETURN
ENDIF
CONDIT = 1.0 / RCOND
ELSEIF(SOLVE) THEN
if(dbg) write(lout,'('' photwo - solve'')')
IF(LEVEL(2).GT.1) THEN
write(lout,'('' Residual: '')')
write(lout,'(1x,1p6e10.3)') (buffer(i),i=1,natj)
ENDIF
C*****precision > double
CALL DGESL (A, NATJ, NATJ, IPVT, BUFFER, 0)
C*****END precision > double
C*****precision > single
C CALL SGESL (A, NATJ, NATJ, IPVT, BUFFER, 0)
C*****END precision > single
IF(LEVEL(2).GT.1) THEN
write(lout,'('' Correction: '')')
write(lout,'(1x,1p6e10.3)') (buffer(i),i=1,natj)
ENDIF
ELSEIF(STORE) THEN
if(dbg) write(lout,'('' photwo - store'')')
C*****precision > double
CALL DCOPY (NATJ, BUFFER, 1, SN, 1)
C*****END precision > double
C*****precision > single
C CALL SCOPY (NATJ, BUFFER, 1, SN, 1)
C*****END precision > single
ELSEIF(SHOW) THEN
write(lout,'('' Solution: '')')
write(lout,'(1x,1p6e10.3)') (s(i),i=1,natj)
ELSEIF(SAVE) THEN
ELSEIF(UPDATE) THEN
CALL WRIT40(' PHOTWO - Should not be trying to update')
LERR=.TRUE.
RETURN
ENDIF
GO TO 500
ENDIF
IF(.NOT. SUCCES .AND. .NOT.LERR) THEN
IF(TPSEM)
1 CALL WRIT40(' PHOTWO - TWOPNT failed to converge ')
CALL SET_ERR(292,
* 'Error. PHOTWO - TWOPNT failed to converge.',1)
C CALL EAROUT(1)
ENDIF
1000 CONTINUE
C
if(dbg) call banner(2,'photwo',0)
C
NAMSUB='photwo'
RETURN
END
C
C------------------------------------------------------------
c
SUBROUTINE PHOJAC (KK, NATJ, LENRGY, ABSOL, RELAT, P, T, VOL,
1 ICKWRK, RCKWRK, WT, SCRTCH, SN, S,
2 RES, RESN, A, SU, AP, H0K, LTIME, DT, T0)
C
C*****precision > double
IMPLICIT DOUBLE PRECISION (A-H, O-Z), INTEGER (I-N)
C*****END precision > double
C*****precision > single
C IMPLICIT REAL (A-H,O-Z), INTEGER (I-N)
C*****END precision > single
C
DIMENSION S(*), RES(*), RESN(*), A(NATJ,*), SN(*), SCRTCH(KK,*),
1 WT(*), ICKWRK(*), RCKWRK(*), SU(*), AP(*), H0K(*)
LOGICAL LENRGY, LTIME
C
COMMON /NAMFN/NAMSUB,NAMFUN
C
CHARACTER*6 NAMSUB,NAMFUN
NAMSUB='PHOJAC'
C
C ZERO THE MATRIX STORAGE SPACE.
C
ZERO = 0.0
C*****precision > double
CALL DCOPY (NATJ * NATJ, ZERO, 0, A, 1)
C*****END precision > double
C
C*****precision > single
C CALL SCOPY (NATJ * NATJ, ZERO, 0, A, 1)
C*****END precision > single
C
C CALL THE FUNCTION AT S AND STORE IN FN.
C
CALL PHOFUN (KK, NATJ, LENRGY, P, T, VOL, ICKWRK, RCKWRK, WT,
2 SCRTCH(1,1), SCRTCH(1,2), SN, S, RESN, SU, AP,
3 H0K, QDOT, LTIME, DT, T0)
C
C TOP OF THE LOOPS OVER THE RESIDUE CLASSES AND
C SOLUTION COMPONENTS.
C
DO 0200 M = 1, NATJ
C
C FOR A GIVEN RESIDUE CLASS AND A GIVEN SOLUTION COMPONENT,
C PERTRB THE S VECTOR AT POINTS IN THE SAME RESIDUE CLASS.
C
SAVE = S(M)
PERTRB = ABS(S(M)) * RELAT + ABSOL
S(M) = S(M) + PERTRB
C
C CALL THE FUNCTION AT THE PERTURBED S AND STORE
C THE RESULT IN F.
C
CALL PHOFUN (KK, NATJ, LENRGY, P, T, VOL, ICKWRK, RCKWRK,
1 WT, SCRTCH(1,1), SCRTCH(1,2), SN, S, RES,
2 SU, AP, H0K, QDOT, LTIME, DT, T0)
C
C RESTORE S TO ITS ORIGINAL VALUE.
C
S(M) = SAVE
C
C DIFFERENCE TO GET THE COLUMNS OF THE JACOBIAN.
C
DO 0100 N = 1, NATJ
A(N, M) = (RES(N) - RESN(N)) / PERTRB
100 CONTINUE
C
C BOTTOM OF THE LOOPS OVER THE RESIDUE CLASSES AND SOLUTION
C COMPONENTS.
C
200 CONTINUE
C
NAMSUB='phojac'
RETURN
END
C
C--------------------------------------------------------------------
c
SUBROUTINE PHOFUN (KK, NATJ, LENRGY, P, T, VOL, ICKWRK, RCKWRK,
1 WT, WDOT, H, SN, S, RES, SU, AP, H0K, QDOT,
2 LTIME, DT, T0)
C
C*****precision > double
IMPLICIT DOUBLE PRECISION (A-H, O-Z), INTEGER (I-N)
C*****END precision > double
C*****precision > single
C IMPLICIT REAL (A-H,O-Z), INTEGER (I-N)
C*****END precision > single
COMMON/CKMCL0/VARPRS,CONTEM
LOGICAL VARPRS,CONTEM
c COMMON/RDATA/TINY,RD1(84)
C
DIMENSION S(*), SN(*), RES(*), WT(*), H(*), WDOT(*),
1 ICKWRK(*), RCKWRK(*), SU(*), AP(*), H0K(*)
C
LOGICAL LENRGY, LTIME
C
CHARACTER*6 NAMSUB,NAMFUN
COMMON /NAMFN/NAMSUB,NAMFUN
DATA TINY/1.D-30/
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C KK - NUMBER OF SPECIES IN THE CHEMICAL SCHEME
C NATJ - NUMBER OF VARIABLES SOLVED
C P - THE PRESSURE.
C CGS UNITS - DYNES/CM**2
C T - THE TEMPERATURE. FOR FIXED-TEMPERATURE PROBLEMS THIS
C IS THE USER-SPECIFIED TEMPERATURE.
C CGS UNITS - K
C VOL - THE VOLUME OF THE CELL
C CGS UNITS - CM**3
C WT - THE ARRAY OF SPECIES MOLECULAR WEIGHTS.
C CGS UNITS - GM/MOLE
C DIMENSION WT(*) AT LEAST KK.
C SN - ARRAY OF SPECIES MASS FRACTIONS PRIOR TO SOLUTION
C DIMENSION SN(*) AT LEAST KK.
C S - DEPENDENT VARIABLE ARRAY. THE TEMPERATURE IS STORED IN
C T=S(NT), AND THE MASS FRACTIONS ARE IN Y(K)=S(NYS+K)
C DIMENSION S(*) AT LEAST NATJ.
C SU - RESIDUALS FOR THE SOLVED VARIABLES AS COMPUTED BY
C PHOENICS (IE. WITHOUT ANY CONTRIBUTION FROM CHEMKIN
C CHEMICAL REACTIONS)
C CGS UNITS - G/SEC (OR ERGS/SEC)
C DIMENSION SU(*) AT LEAST NATJ
C AP - COEFFICIENTS FOR THE SOLVED VARIABLES AS COMPUTED BY
C PHOENICS (IE. WITHOUT ANY CONTRIBUTION FROM CHEMKIN
C CHEMICAL REACTIONS)
C CGS UNITS - G/SEC (OR ERGS/SEC/KELVIN)
C DIMENSION AP(*) AT LEAST NATJ
C H0K - ARRAY OF SPECIES ENTHALPIES AT REFERENCE TEMPERATURE
C CGS UNITS - ERGS/MOLE
C DIMENSION H(*) AT LEAST KK.
C LTIME - IF LTIME=.TRUE. A TIME STEP OF DT WILL BE ADDED INTO
C THE RESIDUAL.
C DT - THE TIME STEP THAT IS USED IF LTIME=.TRUE.
C CGS UNITS - SEC
C
C WORK AND SCRATCH SPACE-
C RCKWRK - FLOATING POINT CHEMKIN WORK SPACE.
C DIMENSIONING - SEE CHEMKIN DOCUMENTATION.
C WDOT - ARRAY OF CHEMICAL PRODUCTION RATES.
C CGS UNITS - MOLES/(CM**3-SEC)
C DIMENSION WDOT(*) AT LEAST KK.
C H - ARRAY OF SPECIES ENTHALPIES.
C CGS UNITS - ERGS/MOLE
C DIMENSION H(*) AT LEAST KK.
C
C OUTPUT-
C RES - RESIDUALS OF THE GOVERNING EQUATIONS AS EVALUATED AT
C S(N). THE RESIDUAL OF THE KK SPECIES EQUATIONS IS IN
C RES(NYS+K), THE ENERGY EQUATION IN RES(NT,J). IF
C INERGY=.FALSE. THEN THE ENERGY EQUATION IS REPLACED
C BY THE GIVEN TEMPERATURE.
C DIMENSION RES(*) AT LEAST NATJ.
C QDOT - HEAT RELEASE RATE FOR THE CURRENT CELL
C CGS UNITS - ERGS/SEC/CC
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C*****unicos timelimit
C call tremain (seconds)
C IF(seconds .lt. 60.0) STOP
C*****END unicos timelimit
NAMSUB='PHOFUN'
C
C FORM THE CHEMICAL RATE TERMS
C
IF(NATJ.EQ.KK) THEN
NY0=1
TL=S(1)+T0
ELSE
NY0=0
IF(CONTEM) THEN
TL=T0
ELSE
TL=T+T0
ENDIF
ENDIF
SUMY=0.0
DO 10 K=1,KK-1
SUMY=SUMY+S(NY0+K)
10 CONTINUE
S(NY0+KK)=1.-SUMY
CALL CKWYP (P, TL, S(NY0+1), ICKWRK, RCKWRK, WDOT)
C
C SPECIES CONSERVATION EQUATION - SU contains the contribution from
C transport terms and sources other than chemical reactions
C
K1=NY0
DO 100 K = 1, KK-1
K1=K1+1
RES(K1) = SU(K1) - AP(K1) * (S(K1) - SN(K1))
1 + VOL * WT(K) * WDOT(K)
100 CONTINUE
IF(LTIME) THEN
CALL CKRHOY (P, TL, S(NY0+1), ICKWRK, RCKWRK, RHO)
K1=NY0
DO 101 K = 1, KK-1
K1=K1+1
RES(K1) = RES(K1) - VOL * RHO * (S(K1) - SN(K1)) / DT
101 CONTINUE
ENDIF
C
C ENERGY EQUATION
C
QDOT=0.0
IF(NATJ.EQ.KK) THEN
IF(LENRGY) THEN
DO 200 K = 1, KK
DQDOT = H0K(K) * WT(K) * WDOT(K)
QDOT = QDOT - DQDOT
200 CONTINUE
ap(1)=ap(1) + tiny
RES(1) = S(1) - SN(1) - (SU(1) + VOL * QDOT) / AP(1)
IF(LTIME) THEN
CALL CKCPBS (TL, S(NY0+1), ICKWRK, RCKWRK, CPMIX)
RES(1) = RES(1) + VOL * RHO * CPMIX * (S(1) - SN(1)) /
1 (DT*AP(1))
ENDIF
ELSE
RES(1)=S(1) - T
ENDIF
ENDIF
C
NAMSUB='phofun'
END
c