|
3911
|
1 C Work performed under the auspices of the U.S. Department of Energy |
|
|
2 C by Lawrence Livermore National Laboratory under contract number |
|
|
3 C W-7405-Eng-48. |
|
|
4 C |
|
|
5 SUBROUTINE DDSTP(X,Y,YPRIME,NEQ,RES,JAC,PSOL,H,WT,VT, |
|
|
6 * JSTART,IDID,RPAR,IPAR,PHI,SAVR,DELTA,E,WM,IWM, |
|
|
7 * ALPHA,BETA,GAMMA,PSI,SIGMA,CJ,CJOLD,HOLD,S,HMIN,UROUND, |
|
|
8 * EPLI,SQRTN,RSQRTN,EPCON,IPHASE,JCALC,JFLG,K,KOLD,NS,NONNEG, |
|
|
9 * NTYPE,NLS) |
|
|
10 C |
|
|
11 C***BEGIN PROLOGUE DDSTP |
|
|
12 C***REFER TO DDASPK |
|
|
13 C***DATE WRITTEN 890101 (YYMMDD) |
|
|
14 C***REVISION DATE 900926 (YYMMDD) |
|
|
15 C***REVISION DATE 940909 (YYMMDD) (Reset PSI(1), PHI(*,2) at 690) |
|
|
16 C |
|
|
17 C |
|
|
18 C----------------------------------------------------------------------- |
|
|
19 C***DESCRIPTION |
|
|
20 C |
|
|
21 C DDSTP solves a system of differential/algebraic equations of |
|
|
22 C the form G(X,Y,YPRIME) = 0, for one step (normally from X to X+H). |
|
|
23 C |
|
|
24 C The methods used are modified divided difference, fixed leading |
|
|
25 C coefficient forms of backward differentiation formulas. |
|
|
26 C The code adjusts the stepsize and order to control the local error |
|
|
27 C per step. |
|
|
28 C |
|
|
29 C |
|
|
30 C The parameters represent |
|
|
31 C X -- Independent variable. |
|
|
32 C Y -- Solution vector at X. |
|
|
33 C YPRIME -- Derivative of solution vector |
|
|
34 C after successful step. |
|
|
35 C NEQ -- Number of equations to be integrated. |
|
|
36 C RES -- External user-supplied subroutine |
|
|
37 C to evaluate the residual. See RES description |
|
|
38 C in DDASPK prologue. |
|
|
39 C JAC -- External user-supplied routine to update |
|
|
40 C Jacobian or preconditioner information in the |
|
|
41 C nonlinear solver. See JAC description in DDASPK |
|
|
42 C prologue. |
|
|
43 C PSOL -- External user-supplied routine to solve |
|
|
44 C a linear system using preconditioning. |
|
|
45 C (This is optional). See PSOL in DDASPK prologue. |
|
|
46 C H -- Appropriate step size for next step. |
|
|
47 C Normally determined by the code. |
|
|
48 C WT -- Vector of weights for error criterion used in Newton test. |
|
|
49 C VT -- Masked vector of weights used in error test. |
|
|
50 C JSTART -- Integer variable set 0 for |
|
|
51 C first step, 1 otherwise. |
|
|
52 C IDID -- Completion code returned from the nonlinear solver. |
|
|
53 C See IDID description in DDASPK prologue. |
|
|
54 C RPAR,IPAR -- Real and integer parameter arrays that |
|
|
55 C are used for communication between the |
|
|
56 C calling program and external user routines. |
|
|
57 C They are not altered by DNSK |
|
|
58 C PHI -- Array of divided differences used by |
|
|
59 C DDSTP. The length is NEQ*(K+1), where |
|
|
60 C K is the maximum order. |
|
|
61 C SAVR -- Work vector for DDSTP of length NEQ. |
|
|
62 C DELTA,E -- Work vectors for DDSTP of length NEQ. |
|
|
63 C WM,IWM -- Real and integer arrays storing |
|
|
64 C information required by the linear solver. |
|
|
65 C |
|
|
66 C The other parameters are information |
|
|
67 C which is needed internally by DDSTP to |
|
|
68 C continue from step to step. |
|
|
69 C |
|
|
70 C----------------------------------------------------------------------- |
|
|
71 C***ROUTINES CALLED |
|
|
72 C NLS, DDWNRM, DDATRP |
|
|
73 C |
|
|
74 C***END PROLOGUE DDSTP |
|
|
75 C |
|
|
76 C |
|
|
77 IMPLICIT DOUBLE PRECISION(A-H,O-Z) |
|
|
78 DIMENSION Y(*),YPRIME(*),WT(*),VT(*) |
|
|
79 DIMENSION PHI(NEQ,*),SAVR(*),DELTA(*),E(*) |
|
|
80 DIMENSION WM(*),IWM(*) |
|
|
81 DIMENSION PSI(*),ALPHA(*),BETA(*),GAMMA(*),SIGMA(*) |
|
|
82 DIMENSION RPAR(*),IPAR(*) |
|
|
83 EXTERNAL RES, JAC, PSOL, NLS |
|
|
84 C |
|
|
85 PARAMETER (LMXORD=3) |
|
|
86 PARAMETER (LNST=11, LETF=14, LCFN=15) |
|
|
87 C |
|
|
88 C |
|
|
89 C----------------------------------------------------------------------- |
|
|
90 C BLOCK 1. |
|
|
91 C Initialize. On the first call, set |
|
|
92 C the order to 1 and initialize |
|
|
93 C other variables. |
|
|
94 C----------------------------------------------------------------------- |
|
|
95 C |
|
|
96 C Initializations for all calls |
|
|
97 C |
|
|
98 XOLD=X |
|
|
99 NCF=0 |
|
|
100 NEF=0 |
|
|
101 IF(JSTART .NE. 0) GO TO 120 |
|
|
102 C |
|
|
103 C If this is the first step, perform |
|
|
104 C other initializations |
|
|
105 C |
|
|
106 K=1 |
|
|
107 KOLD=0 |
|
|
108 HOLD=0.0D0 |
|
|
109 PSI(1)=H |
|
|
110 CJ = 1.D0/H |
|
|
111 IPHASE = 0 |
|
|
112 NS=0 |
|
|
113 120 CONTINUE |
|
|
114 C |
|
|
115 C |
|
|
116 C |
|
|
117 C |
|
|
118 C |
|
|
119 C----------------------------------------------------------------------- |
|
|
120 C BLOCK 2 |
|
|
121 C Compute coefficients of formulas for |
|
|
122 C this step. |
|
|
123 C----------------------------------------------------------------------- |
|
|
124 200 CONTINUE |
|
|
125 KP1=K+1 |
|
|
126 KP2=K+2 |
|
|
127 KM1=K-1 |
|
|
128 IF(H.NE.HOLD.OR.K .NE. KOLD) NS = 0 |
|
|
129 NS=MIN0(NS+1,KOLD+2) |
|
|
130 NSP1=NS+1 |
|
|
131 IF(KP1 .LT. NS)GO TO 230 |
|
|
132 C |
|
|
133 BETA(1)=1.0D0 |
|
|
134 ALPHA(1)=1.0D0 |
|
|
135 TEMP1=H |
|
|
136 GAMMA(1)=0.0D0 |
|
|
137 SIGMA(1)=1.0D0 |
|
|
138 DO 210 I=2,KP1 |
|
|
139 TEMP2=PSI(I-1) |
|
|
140 PSI(I-1)=TEMP1 |
|
|
141 BETA(I)=BETA(I-1)*PSI(I-1)/TEMP2 |
|
|
142 TEMP1=TEMP2+H |
|
|
143 ALPHA(I)=H/TEMP1 |
|
|
144 SIGMA(I)=(I-1)*SIGMA(I-1)*ALPHA(I) |
|
|
145 GAMMA(I)=GAMMA(I-1)+ALPHA(I-1)/H |
|
|
146 210 CONTINUE |
|
|
147 PSI(KP1)=TEMP1 |
|
|
148 230 CONTINUE |
|
|
149 C |
|
|
150 C Compute ALPHAS, ALPHA0 |
|
|
151 C |
|
|
152 ALPHAS = 0.0D0 |
|
|
153 ALPHA0 = 0.0D0 |
|
|
154 DO 240 I = 1,K |
|
|
155 ALPHAS = ALPHAS - 1.0D0/I |
|
|
156 ALPHA0 = ALPHA0 - ALPHA(I) |
|
|
157 240 CONTINUE |
|
|
158 C |
|
|
159 C Compute leading coefficient CJ |
|
|
160 C |
|
|
161 CJLAST = CJ |
|
|
162 CJ = -ALPHAS/H |
|
|
163 C |
|
|
164 C Compute variable stepsize error coefficient CK |
|
|
165 C |
|
|
166 CK = ABS(ALPHA(KP1) + ALPHAS - ALPHA0) |
|
|
167 CK = MAX(CK,ALPHA(KP1)) |
|
|
168 C |
|
|
169 C Change PHI to PHI STAR |
|
|
170 C |
|
|
171 IF(KP1 .LT. NSP1) GO TO 280 |
|
|
172 DO 270 J=NSP1,KP1 |
|
|
173 DO 260 I=1,NEQ |
|
|
174 260 PHI(I,J)=BETA(J)*PHI(I,J) |
|
|
175 270 CONTINUE |
|
|
176 280 CONTINUE |
|
|
177 C |
|
|
178 C Update time |
|
|
179 C |
|
|
180 X=X+H |
|
|
181 C |
|
|
182 C Initialize IDID to 1 |
|
|
183 C |
|
|
184 IDID = 1 |
|
|
185 C |
|
|
186 C |
|
|
187 C |
|
|
188 C |
|
|
189 C |
|
|
190 C----------------------------------------------------------------------- |
|
|
191 C BLOCK 3 |
|
|
192 C Call the nonlinear system solver to obtain the solution and |
|
|
193 C derivative. |
|
|
194 C----------------------------------------------------------------------- |
|
|
195 C |
|
|
196 CALL NLS(X,Y,YPRIME,NEQ, |
|
|
197 * RES,JAC,PSOL,H,WT,JSTART,IDID,RPAR,IPAR,PHI,GAMMA, |
|
|
198 * SAVR,DELTA,E,WM,IWM,CJ,CJOLD,CJLAST,S, |
|
|
199 * UROUND,EPLI,SQRTN,RSQRTN,EPCON,JCALC,JFLG,KP1, |
|
|
200 * NONNEG,NTYPE,IERNLS) |
|
|
201 C |
|
|
202 IF(IERNLS .NE. 0)GO TO 600 |
|
|
203 C |
|
|
204 C |
|
|
205 C |
|
|
206 C |
|
|
207 C |
|
|
208 C----------------------------------------------------------------------- |
|
|
209 C BLOCK 4 |
|
|
210 C Estimate the errors at orders K,K-1,K-2 |
|
|
211 C as if constant stepsize was used. Estimate |
|
|
212 C the local error at order K and test |
|
|
213 C whether the current step is successful. |
|
|
214 C----------------------------------------------------------------------- |
|
|
215 C |
|
|
216 C Estimate errors at orders K,K-1,K-2 |
|
|
217 C |
|
|
218 ENORM = DDWNRM(NEQ,E,VT,RPAR,IPAR) |
|
|
219 ERK = SIGMA(K+1)*ENORM |
|
|
220 TERK = (K+1)*ERK |
|
|
221 EST = ERK |
|
|
222 KNEW=K |
|
|
223 IF(K .EQ. 1)GO TO 430 |
|
|
224 DO 405 I = 1,NEQ |
|
|
225 405 DELTA(I) = PHI(I,KP1) + E(I) |
|
|
226 ERKM1=SIGMA(K)*DDWNRM(NEQ,DELTA,VT,RPAR,IPAR) |
|
|
227 TERKM1 = K*ERKM1 |
|
|
228 IF(K .GT. 2)GO TO 410 |
|
|
229 IF(TERKM1 .LE. 0.5*TERK)GO TO 420 |
|
|
230 GO TO 430 |
|
|
231 410 CONTINUE |
|
|
232 DO 415 I = 1,NEQ |
|
|
233 415 DELTA(I) = PHI(I,K) + DELTA(I) |
|
|
234 ERKM2=SIGMA(K-1)*DDWNRM(NEQ,DELTA,VT,RPAR,IPAR) |
|
|
235 TERKM2 = (K-1)*ERKM2 |
|
|
236 IF(MAX(TERKM1,TERKM2).GT.TERK)GO TO 430 |
|
|
237 C |
|
|
238 C Lower the order |
|
|
239 C |
|
|
240 420 CONTINUE |
|
|
241 KNEW=K-1 |
|
|
242 EST = ERKM1 |
|
|
243 C |
|
|
244 C |
|
|
245 C Calculate the local error for the current step |
|
|
246 C to see if the step was successful |
|
|
247 C |
|
|
248 430 CONTINUE |
|
|
249 ERR = CK * ENORM |
|
|
250 IF(ERR .GT. 1.0D0)GO TO 600 |
|
|
251 C |
|
|
252 C |
|
|
253 C |
|
|
254 C |
|
|
255 C |
|
|
256 C----------------------------------------------------------------------- |
|
|
257 C BLOCK 5 |
|
|
258 C The step is successful. Determine |
|
|
259 C the best order and stepsize for |
|
|
260 C the next step. Update the differences |
|
|
261 C for the next step. |
|
|
262 C----------------------------------------------------------------------- |
|
|
263 IDID=1 |
|
|
264 IWM(LNST)=IWM(LNST)+1 |
|
|
265 KDIFF=K-KOLD |
|
|
266 KOLD=K |
|
|
267 HOLD=H |
|
|
268 C |
|
|
269 C |
|
|
270 C Estimate the error at order K+1 unless |
|
|
271 C already decided to lower order, or |
|
|
272 C already using maximum order, or |
|
|
273 C stepsize not constant, or |
|
|
274 C order raised in previous step |
|
|
275 C |
|
|
276 IF(KNEW.EQ.KM1.OR.K.EQ.IWM(LMXORD))IPHASE=1 |
|
|
277 IF(IPHASE .EQ. 0)GO TO 545 |
|
|
278 IF(KNEW.EQ.KM1)GO TO 540 |
|
|
279 IF(K.EQ.IWM(LMXORD)) GO TO 550 |
|
|
280 IF(KP1.GE.NS.OR.KDIFF.EQ.1)GO TO 550 |
|
|
281 DO 510 I=1,NEQ |
|
|
282 510 DELTA(I)=E(I)-PHI(I,KP2) |
|
|
283 ERKP1 = (1.0D0/(K+2))*DDWNRM(NEQ,DELTA,VT,RPAR,IPAR) |
|
|
284 TERKP1 = (K+2)*ERKP1 |
|
|
285 IF(K.GT.1)GO TO 520 |
|
|
286 IF(TERKP1.GE.0.5D0*TERK)GO TO 550 |
|
|
287 GO TO 530 |
|
|
288 520 IF(TERKM1.LE.MIN(TERK,TERKP1))GO TO 540 |
|
|
289 IF(TERKP1.GE.TERK.OR.K.EQ.IWM(LMXORD))GO TO 550 |
|
|
290 C |
|
|
291 C Raise order |
|
|
292 C |
|
|
293 530 K=KP1 |
|
|
294 EST = ERKP1 |
|
|
295 GO TO 550 |
|
|
296 C |
|
|
297 C Lower order |
|
|
298 C |
|
|
299 540 K=KM1 |
|
|
300 EST = ERKM1 |
|
|
301 GO TO 550 |
|
|
302 C |
|
|
303 C If IPHASE = 0, increase order by one and multiply stepsize by |
|
|
304 C factor two |
|
|
305 C |
|
|
306 545 K = KP1 |
|
|
307 HNEW = H*2.0D0 |
|
|
308 H = HNEW |
|
|
309 GO TO 575 |
|
|
310 C |
|
|
311 C |
|
|
312 C Determine the appropriate stepsize for |
|
|
313 C the next step. |
|
|
314 C |
|
|
315 550 HNEW=H |
|
|
316 TEMP2=K+1 |
|
|
317 R=(2.0D0*EST+0.0001D0)**(-1.0D0/TEMP2) |
|
|
318 IF(R .LT. 2.0D0) GO TO 555 |
|
|
319 HNEW = 2.0D0*H |
|
|
320 GO TO 560 |
|
|
321 555 IF(R .GT. 1.0D0) GO TO 560 |
|
|
322 R = MAX(0.5D0,MIN(0.9D0,R)) |
|
|
323 HNEW = H*R |
|
|
324 560 H=HNEW |
|
|
325 C |
|
|
326 C |
|
|
327 C Update differences for next step |
|
|
328 C |
|
|
329 575 CONTINUE |
|
|
330 IF(KOLD.EQ.IWM(LMXORD))GO TO 585 |
|
|
331 DO 580 I=1,NEQ |
|
|
332 580 PHI(I,KP2)=E(I) |
|
|
333 585 CONTINUE |
|
|
334 DO 590 I=1,NEQ |
|
|
335 590 PHI(I,KP1)=PHI(I,KP1)+E(I) |
|
|
336 DO 595 J1=2,KP1 |
|
|
337 J=KP1-J1+1 |
|
|
338 DO 595 I=1,NEQ |
|
|
339 595 PHI(I,J)=PHI(I,J)+PHI(I,J+1) |
|
|
340 JSTART = 1 |
|
|
341 RETURN |
|
|
342 C |
|
|
343 C |
|
|
344 C |
|
|
345 C |
|
|
346 C |
|
|
347 C----------------------------------------------------------------------- |
|
|
348 C BLOCK 6 |
|
|
349 C The step is unsuccessful. Restore X,PSI,PHI |
|
|
350 C Determine appropriate stepsize for |
|
|
351 C continuing the integration, or exit with |
|
|
352 C an error flag if there have been many |
|
|
353 C failures. |
|
|
354 C----------------------------------------------------------------------- |
|
|
355 600 IPHASE = 1 |
|
|
356 C |
|
|
357 C Restore X,PHI,PSI |
|
|
358 C |
|
|
359 X=XOLD |
|
|
360 IF(KP1.LT.NSP1)GO TO 630 |
|
|
361 DO 620 J=NSP1,KP1 |
|
|
362 TEMP1=1.0D0/BETA(J) |
|
|
363 DO 610 I=1,NEQ |
|
|
364 610 PHI(I,J)=TEMP1*PHI(I,J) |
|
|
365 620 CONTINUE |
|
|
366 630 CONTINUE |
|
|
367 DO 640 I=2,KP1 |
|
|
368 640 PSI(I-1)=PSI(I)-H |
|
|
369 C |
|
|
370 C |
|
|
371 C Test whether failure is due to nonlinear solver |
|
|
372 C or error test |
|
|
373 C |
|
|
374 IF(IERNLS .EQ. 0)GO TO 660 |
|
|
375 IWM(LCFN)=IWM(LCFN)+1 |
|
|
376 C |
|
|
377 C |
|
|
378 C The nonlinear solver failed to converge. |
|
|
379 C Determine the cause of the failure and take appropriate action. |
|
|
380 C If IERNLS .LT. 0, then return. Otherwise, reduce the stepsize |
|
|
381 C and try again, unless too many failures have occurred. |
|
|
382 C |
|
|
383 IF (IERNLS .LT. 0) GO TO 675 |
|
|
384 NCF = NCF + 1 |
|
|
385 R = 0.25D0 |
|
|
386 H = H*R |
|
|
387 IF (NCF .LT. 10 .AND. ABS(H) .GE. HMIN) GO TO 690 |
|
|
388 IF (IDID .EQ. 1) IDID = -7 |
|
|
389 IF (NEF .GE. 3) IDID = -9 |
|
|
390 GO TO 675 |
|
|
391 C |
|
|
392 C |
|
|
393 C The nonlinear solver converged, and the cause |
|
|
394 C of the failure was the error estimate |
|
|
395 C exceeding the tolerance. |
|
|
396 C |
|
|
397 660 NEF=NEF+1 |
|
|
398 IWM(LETF)=IWM(LETF)+1 |
|
|
399 IF (NEF .GT. 1) GO TO 665 |
|
|
400 C |
|
|
401 C On first error test failure, keep current order or lower |
|
|
402 C order by one. Compute new stepsize based on differences |
|
|
403 C of the solution. |
|
|
404 C |
|
|
405 K = KNEW |
|
|
406 TEMP2 = K + 1 |
|
|
407 R = 0.90D0*(2.0D0*EST+0.0001D0)**(-1.0D0/TEMP2) |
|
|
408 R = MAX(0.25D0,MIN(0.9D0,R)) |
|
|
409 H = H*R |
|
|
410 IF (ABS(H) .GE. HMIN) GO TO 690 |
|
|
411 IDID = -6 |
|
|
412 GO TO 675 |
|
|
413 C |
|
|
414 C On second error test failure, use the current order or |
|
|
415 C decrease order by one. Reduce the stepsize by a factor of |
|
|
416 C one quarter. |
|
|
417 C |
|
|
418 665 IF (NEF .GT. 2) GO TO 670 |
|
|
419 K = KNEW |
|
|
420 R = 0.25D0 |
|
|
421 H = R*H |
|
|
422 IF (ABS(H) .GE. HMIN) GO TO 690 |
|
|
423 IDID = -6 |
|
|
424 GO TO 675 |
|
|
425 C |
|
|
426 C On third and subsequent error test failures, set the order to |
|
|
427 C one, and reduce the stepsize by a factor of one quarter. |
|
|
428 C |
|
|
429 670 K = 1 |
|
|
430 R = 0.25D0 |
|
|
431 H = R*H |
|
|
432 IF (ABS(H) .GE. HMIN) GO TO 690 |
|
|
433 IDID = -6 |
|
|
434 GO TO 675 |
|
|
435 C |
|
|
436 C |
|
|
437 C |
|
|
438 C |
|
|
439 C For all crashes, restore Y to its last value, |
|
|
440 C interpolate to find YPRIME at last X, and return. |
|
|
441 C |
|
|
442 C Before returning, verify that the user has not set |
|
|
443 C IDID to a nonnegative value. If the user has set IDID |
|
|
444 C to a nonnegative value, then reset IDID to be -7, indicating |
|
|
445 C a failure in the nonlinear system solver. |
|
|
446 C |
|
|
447 675 CONTINUE |
|
|
448 CALL DDATRP(X,X,Y,YPRIME,NEQ,K,PHI,PSI) |
|
|
449 JSTART = 1 |
|
|
450 IF (IDID .GE. 0) IDID = -7 |
|
|
451 RETURN |
|
|
452 C |
|
|
453 C |
|
|
454 C Go back and try this step again. |
|
|
455 C If this is the first step, reset PSI(1) and rescale PHI(*,2). |
|
|
456 C |
|
|
457 690 IF (KOLD .EQ. 0) THEN |
|
|
458 PSI(1) = H |
|
|
459 DO 695 I = 1,NEQ |
|
|
460 695 PHI(I,2) = R*PHI(I,2) |
|
|
461 ENDIF |
|
|
462 GO TO 200 |
|
|
463 C |
|
|
464 C------END OF SUBROUTINE DDSTP------------------------------------------ |
|
|
465 END |
