cgxbfgr.for

c
c

File-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) 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',190495) call writ2r('ap ',ap,'su ',su) 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