00001 C Work performed under the auspices of the U.S. Department of Energy 00002 C by Lawrence Livermore National Laboratory under contract number 00003 C W-7405-Eng-48. 00004 C 00005 SUBROUTINE DNEDK(X,Y,YPRIME,NEQ,RES,JACK,PSOL, 00006 * H,WT,JSTART,IDID,RPAR,IPAR,PHI,GAMMA,SAVR,DELTA,E, 00007 * WM,IWM,CJ,CJOLD,CJLAST,S,UROUND,EPLI,SQRTN,RSQRTN, 00008 * EPCON,JCALC,JFLG,KP1,NONNEG,NTYPE,IERNLS) 00009 C 00010 C***BEGIN PROLOGUE DNEDK 00011 C***REFER TO DDASPK 00012 C***DATE WRITTEN 891219 (YYMMDD) 00013 C***REVISION DATE 900926 (YYMMDD) 00014 C***REVISION DATE 940701 (YYMMDD) 00015 C 00016 C 00017 C----------------------------------------------------------------------- 00018 C***DESCRIPTION 00019 C 00020 C DNEDK solves a nonlinear system of 00021 C algebraic equations of the form 00022 C G(X,Y,YPRIME) = 0 for the unknown Y. 00023 C 00024 C The method used is a matrix-free Newton scheme. 00025 C 00026 C The parameters represent 00027 C X -- Independent variable. 00028 C Y -- Solution vector at x. 00029 C YPRIME -- Derivative of solution vector 00030 C after successful step. 00031 C NEQ -- Number of equations to be integrated. 00032 C RES -- External user-supplied subroutine 00033 C to evaluate the residual. See RES description 00034 C in DDASPK prologue. 00035 C JACK -- External user-supplied routine to update 00036 C the preconditioner. (This is optional). 00037 C See JAC description for the case 00038 C INFO(12) = 1 in the DDASPK prologue. 00039 C PSOL -- External user-supplied routine to solve 00040 C a linear system using preconditioning. 00041 C (This is optional). See explanation inside DDASPK. 00042 C H -- Appropriate step size for this step. 00043 C WT -- Vector of weights for error criterion. 00044 C JSTART -- Indicates first call to this routine. 00045 C If JSTART = 0, then this is the first call, 00046 C otherwise it is not. 00047 C IDID -- Completion flag, output by DNEDK. 00048 C See IDID description in DDASPK prologue. 00049 C RPAR,IPAR -- Real and integer arrays used for communication 00050 C between the calling program and external user 00051 C routines. They are not altered within DASPK. 00052 C PHI -- Array of divided differences used by 00053 C DNEDK. The length is NEQ*(K+1), where 00054 C K is the maximum order. 00055 C GAMMA -- Array used to predict Y and YPRIME. The length 00056 C is K+1, where K is the maximum order. 00057 C SAVR -- Work vector for DNEDK of length NEQ. 00058 C DELTA -- Work vector for DNEDK of length NEQ. 00059 C E -- Error accumulation vector for DNEDK of length NEQ. 00060 C WM,IWM -- Real and integer arrays storing 00061 C matrix information for linear system 00062 C solvers, and various other information. 00063 C CJ -- Parameter always proportional to 1/H. 00064 C CJOLD -- Saves the value of CJ as of the last call to DITMD. 00065 C Accounts for changes in CJ needed to 00066 C decide whether to call DITMD. 00067 C CJLAST -- Previous value of CJ. 00068 C S -- A scalar determined by the approximate rate 00069 C of convergence of the Newton iteration and used 00070 C in the convergence test for the Newton iteration. 00071 C 00072 C If RATE is defined to be an estimate of the 00073 C rate of convergence of the Newton iteration, 00074 C then S = RATE/(1.D0-RATE). 00075 C 00076 C The closer RATE is to 0., the faster the Newton 00077 C iteration is converging; the closer RATE is to 1., 00078 C the slower the Newton iteration is converging. 00079 C 00080 C On the first Newton iteration with an up-dated 00081 C preconditioner S = 100.D0, Thus the initial 00082 C RATE of convergence is approximately 1. 00083 C 00084 C S is preserved from call to call so that the rate 00085 C estimate from a previous step can be applied to 00086 C the current step. 00087 C UROUND -- Unit roundoff. 00088 C EPLI -- convergence test constant. 00089 C See DDASPK prologue for more details. 00090 C SQRTN -- Square root of NEQ. 00091 C RSQRTN -- reciprical of square root of NEQ. 00092 C EPCON -- Tolerance to test for convergence of the Newton 00093 C iteration. 00094 C JCALC -- Flag used to determine when to update 00095 C the Jacobian matrix. In general: 00096 C 00097 C JCALC = -1 ==> Call the DITMD routine to update 00098 C the Jacobian matrix. 00099 C JCALC = 0 ==> Jacobian matrix is up-to-date. 00100 C JCALC = 1 ==> Jacobian matrix is out-dated, 00101 C but DITMD will not be called unless 00102 C JCALC is set to -1. 00103 C JFLG -- Flag showing whether a Jacobian routine is supplied. 00104 C KP1 -- The current order + 1; updated across calls. 00105 C NONNEG -- Flag to determine nonnegativity constraints. 00106 C NTYPE -- Identification code for the DNEDK routine. 00107 C 1 ==> modified Newton; iterative linear solver. 00108 C 2 ==> modified Newton; user-supplied linear solver. 00109 C IERNLS -- Error flag for nonlinear solver. 00110 C 0 ==> nonlinear solver converged. 00111 C 1 ==> recoverable error inside non-linear solver. 00112 C -1 ==> unrecoverable error inside non-linear solver. 00113 C 00114 C The following group of variables are passed as arguments to 00115 C the Newton iteration solver. They are explained in greater detail 00116 C in DNSK: 00117 C TOLNEW, MULDEL, MAXIT, IERNEW 00118 C 00119 C IERTYP -- Flag which tells whether this subroutine is correct. 00120 C 0 ==> correct subroutine. 00121 C 1 ==> incorrect subroutine. 00122 C 00123 C----------------------------------------------------------------------- 00124 C***ROUTINES CALLED 00125 C RES, JACK, DDWNRM, DNSK 00126 C 00127 C***END PROLOGUE DNEDK 00128 C 00129 C 00130 IMPLICIT DOUBLE PRECISION(A-H,O-Z) 00131 DIMENSION Y(*),YPRIME(*),WT(*) 00132 DIMENSION PHI(NEQ,*),SAVR(*),DELTA(*),E(*) 00133 DIMENSION WM(*),IWM(*) 00134 DIMENSION GAMMA(*),RPAR(*),IPAR(*) 00135 EXTERNAL RES, JACK, PSOL 00136 C 00137 PARAMETER (LNRE=12, LNJE=13, LLOCWP=29, LLCIWP=30) 00138 C 00139 SAVE MULDEL, MAXIT, XRATE 00140 DATA MULDEL/0/, MAXIT/4/, XRATE/0.25D0/ 00141 C 00142 C Verify that this is the correct subroutine. 00143 C 00144 IERTYP = 0 00145 IF (NTYPE .NE. 1) THEN 00146 IERTYP = 1 00147 GO TO 380 00148 ENDIF 00149 C 00150 C If this is the first step, perform initializations. 00151 C 00152 IF (JSTART .EQ. 0) THEN 00153 CJOLD = CJ 00154 JCALC = -1 00155 S = 100.D0 00156 ENDIF 00157 C 00158 C Perform all other initializations. 00159 C 00160 IERNLS = 0 00161 LWP = IWM(LLOCWP) 00162 LIWP = IWM(LLCIWP) 00163 C 00164 C Decide whether to update the preconditioner. 00165 C 00166 IF (JFLG .NE. 0) THEN 00167 TEMP1 = (1.0D0 - XRATE)/(1.0D0 + XRATE) 00168 TEMP2 = 1.0D0/TEMP1 00169 IF (CJ/CJOLD .LT. TEMP1 .OR. CJ/CJOLD .GT. TEMP2) JCALC = -1 00170 IF (CJ .NE. CJLAST) S = 100.D0 00171 ELSE 00172 JCALC = 0 00173 ENDIF 00174 C 00175 C Looping point for updating preconditioner with current stepsize. 00176 C 00177 300 CONTINUE 00178 C 00179 C Initialize all error flags to zero. 00180 C 00181 IERPJ = 0 00182 IRES = 0 00183 IERSL = 0 00184 IERNEW = 0 00185 C 00186 C Predict the solution and derivative and compute the tolerance 00187 C for the Newton iteration. 00188 C 00189 DO 310 I=1,NEQ 00190 Y(I)=PHI(I,1) 00191 310 YPRIME(I)=0.0D0 00192 DO 330 J=2,KP1 00193 DO 320 I=1,NEQ 00194 Y(I)=Y(I)+PHI(I,J) 00195 320 YPRIME(I)=YPRIME(I)+GAMMA(J)*PHI(I,J) 00196 330 CONTINUE 00197 EPLIN = EPLI*EPCON 00198 TOLNEW = EPLIN 00199 C 00200 C Call RES to initialize DELTA. 00201 C 00202 IWM(LNRE)=IWM(LNRE)+1 00203 CALL RES(X,Y,YPRIME,CJ,DELTA,IRES,RPAR,IPAR) 00204 IF (IRES .LT. 0) GO TO 380 00205 C 00206 C 00207 C If indicated, update the preconditioner. 00208 C Set JCALC to 0 as an indicator that this has been done. 00209 C 00210 IF(JCALC .EQ. -1)THEN 00211 IWM(LNJE) = IWM(LNJE) + 1 00212 JCALC=0 00213 CALL JACK (RES, IRES, NEQ, X, Y, YPRIME, WT, DELTA, E, H, CJ, 00214 * WM(LWP), IWM(LIWP), IERPJ, RPAR, IPAR) 00215 CJOLD=CJ 00216 S = 100.D0 00217 IF (IRES .LT. 0) GO TO 380 00218 IF (IERPJ .NE. 0) GO TO 380 00219 ENDIF 00220 C 00221 C Call the nonlinear Newton solver. 00222 C 00223 CALL DNSK(X,Y,YPRIME,NEQ,RES,PSOL,WT,RPAR,IPAR,SAVR, 00224 * DELTA,E,WM,IWM,CJ,SQRTN,RSQRTN,EPLIN,EPCON, 00225 * S,TEMP1,TOLNEW,MULDEL,MAXIT,IRES,IERSL,IERNEW) 00226 C 00227 IF (IERNEW .GT. 0 .AND. JCALC .NE. 0) THEN 00228 C 00229 C The Newton iteration had a recoverable failure with an old 00230 C preconditioner. Retry the step with a new preconditioner. 00231 C 00232 JCALC = -1 00233 GO TO 300 00234 ENDIF 00235 C 00236 IF (IERNEW .NE. 0) GO TO 380 00237 C 00238 C The Newton iteration has converged. If nonnegativity of 00239 C solution is required, set the solution nonnegative, if the 00240 C perturbation to do it is small enough. If the change is too 00241 C large, then consider the corrector iteration to have failed. 00242 C 00243 IF(NONNEG .EQ. 0) GO TO 390 00244 DO 360 I = 1,NEQ 00245 360 DELTA(I) = MIN(Y(I),0.0D0) 00246 DELNRM = DDWNRM(NEQ,DELTA,WT,RPAR,IPAR) 00247 IF(DELNRM .GT. EPCON) GO TO 380 00248 DO 370 I = 1,NEQ 00249 370 E(I) = E(I) - DELTA(I) 00250 GO TO 390 00251 C 00252 C 00253 C Exits from nonlinear solver. 00254 C No convergence with current preconditioner. 00255 C Compute IERNLS and IDID accordingly. 00256 C 00257 380 CONTINUE 00258 IF (IRES .LE. -2 .OR. IERSL .LT. 0 .OR. IERTYP .NE. 0) THEN 00259 IERNLS = -1 00260 IF (IRES .LE. -2) IDID = -11 00261 IF (IERSL .LT. 0) IDID = -13 00262 IF (IERTYP .NE. 0) IDID = -15 00263 ELSE 00264 IERNLS = 1 00265 IF (IRES .EQ. -1) IDID = -10 00266 IF (IERPJ .NE. 0) IDID = -5 00267 IF (IERSL .GT. 0) IDID = -14 00268 ENDIF 00269 C 00270 C 00271 390 JCALC = 1 00272 RETURN 00273 C 00274 C------END OF SUBROUTINE DNEDK------------------------------------------ 00275 END
1.7.1
