Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2019 The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 : module mod_dop853
21 : use const_def, only: dp
22 : use math_lib
23 :
24 : contains
25 :
26 1 : subroutine do_dop853(
27 : & n,fcn,x,y,xend,h,max_step_size,max_steps,
28 1 : & rtol,atol,itol,solout,iout,work,lwork,iwork,liwork,
29 : & lrpar,rpar,lipar,ipar,lout,idid)
30 : ! *** *** *** *** *** *** *** *** *** *** *** *** ***
31 : ! declarations
32 : ! *** *** *** *** *** *** *** *** *** *** *** *** ***
33 : implicit real(dp) (a-h,o-z)
34 : integer, intent(in) :: n ! the dimension of the system
35 : interface
36 : #include "num_fcn.dek"
37 : #include "num_solout.dek"
38 : end interface
39 : real(dp), intent(inout) :: x
40 : real(dp), intent(inout), pointer :: y(:) ! (n)
41 : real(dp), intent(in) :: xend, h, max_step_size
42 : real(dp), intent(in) :: rtol(*)
43 : real(dp), intent(in) :: atol(*)
44 : integer, intent(in) :: itol, iout, liwork, lwork, max_steps
45 : integer, intent(inout) :: iwork(liwork) ! should be 0 on entry
46 : real(dp), intent(inout) :: work(lwork)
47 : integer, intent(in) :: lrpar, lipar
48 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
49 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
50 : integer, intent(in) :: lout
51 : integer, intent(out) :: idid
52 :
53 : logical :: arret
54 : ! *** *** *** *** *** *** ***
55 : ! setting the parameters
56 : ! *** *** *** *** *** *** ***
57 1 : nfcn=0
58 1 : nstep=0
59 1 : naccpt=0
60 1 : nrejct=0
61 1 : arret=.false.
62 : ! -------- nmax , the maximal number of steps -----
63 1 : if(max_steps == 0)then
64 1 : nmax=100000
65 : else
66 0 : nmax=max_steps
67 0 : if(nmax <= 0)then
68 0 : if (lout >0) write(lout,*)
69 0 : & ' wrong input max_steps=',max_steps
70 : arret=.true.
71 : end if
72 : end if
73 : ! -------- meth coefficients of the method
74 1 : if(iwork(2) == 0)then
75 1 : meth=1
76 : else
77 0 : meth=iwork(2)
78 0 : if(meth <= 0.or.meth >= 4)then
79 0 : if (lout >0) write(lout,*)
80 0 : & ' curious input iwork(2)=',iwork(2)
81 : arret=.true.
82 : end if
83 : end if
84 : ! -------- nstiff parameter for stiffness detection
85 1 : nstiff=iwork(4)
86 1 : if (nstiff == 0) nstiff=1000
87 1 : if (nstiff <0) nstiff=nmax+10
88 : ! -------- nrdens number of dense output components
89 1 : nrdens=iwork(5)
90 1 : if(nrdens <0.or.nrdens >n)then
91 0 : if (lout >0) write(lout,*)
92 0 : & ' curious input iwork(5)=',iwork(5)
93 : arret=.true.
94 : else
95 1 : if(nrdens >0.and.iout <2)then
96 0 : if (lout >0) write(lout,*)
97 0 : & ' warning: put iout=2 for dense output '
98 : end if
99 1 : if (nrdens == n) then
100 3 : do i=1,nrdens
101 3 : iwork(i+20)=i
102 : end do
103 : end if
104 : end if
105 : ! -------- uround smallest number satisfying 1.d0+uround>1.d0
106 1 : if(work(1) == 0.d0)then
107 1 : uround=2.3d-16
108 : else
109 0 : uround=work(1)
110 0 : if(uround <= 1.d-35.or.uround >= 1.d0)then
111 0 : if (lout >0) write(lout,*)
112 0 : & ' which machine do you have? your uround was:',work(1)
113 : arret=.true.
114 : end if
115 : end if
116 : ! ------- safety factor -------------
117 1 : if(work(2) == 0.d0)then
118 1 : safe=0.9d0
119 : else
120 0 : safe=work(2)
121 0 : if(safe >= 1.d0.or.safe <= 1.d-4)then
122 0 : if (lout >0) write(lout,*)
123 0 : & ' curious input for safety factor work(2)=',work(2)
124 : arret=.true.
125 : end if
126 : end if
127 : ! ------- fac1,fac2 parameters for step size selection
128 1 : if(work(3) == 0.d0)then
129 1 : fac1=0.333d0
130 : else
131 0 : fac1=work(3)
132 : end if
133 1 : if(work(4) == 0.d0)then
134 1 : fac2=6.d0
135 : else
136 0 : fac2=work(4)
137 : end if
138 : ! --------- beta for step control stabilization -----------
139 1 : if(work(5) == 0.d0)then
140 1 : beta=0.0d0
141 : else
142 0 : if(work(5) <0.d0)then
143 0 : beta=0.d0
144 : else
145 0 : beta=work(5)
146 0 : if(beta >0.2d0)then
147 0 : if (lout >0) write(lout,*)
148 0 : & ' curious input for beta: work(5)=',work(5)
149 : arret=.true.
150 : end if
151 : end if
152 : end if
153 : ! -------- maximal step size
154 1 : if(max_step_size == 0.d0)then
155 1 : hmax=xend-x
156 : else
157 1 : hmax=max_step_size
158 : end if
159 : ! ------- prepare the entry-points for the arrays in work -----
160 1 : iek1=21
161 1 : iek2=iek1+n
162 1 : iek3=iek2+n
163 1 : iek4=iek3+n
164 1 : iek5=iek4+n
165 1 : iek6=iek5+n
166 1 : iek7=iek6+n
167 1 : iek8=iek7+n
168 1 : iek9=iek8+n
169 1 : iek10=iek9+n
170 1 : iey1=iek10+n
171 1 : ieco=iey1+n
172 : ! ------ total storage requirement -----------
173 1 : istore=ieco+(3+8*nrdens)-1
174 1 : if(istore >lwork)then
175 0 : if (lout >0) write(lout,*)
176 0 : & ' insufficient storage for work, min. lwork=',istore
177 : arret=.true.
178 : end if
179 1 : icomp=21
180 1 : istore=icomp+nrdens-1
181 1 : if(istore >liwork)then
182 0 : if (lout >0) write(lout,*)
183 0 : & ' insufficient storage for iwork, min. liwork=',istore
184 : arret=.true.
185 : end if
186 : ! -------- when a fail has occurred, we return with idid=-1
187 1 : if (arret) then
188 0 : idid=-1
189 0 : return
190 : end if
191 : ! -------- call to core integrator ------------
192 : call dp86co(n,fcn,x,y,xend,hmax,h,rtol,atol,itol,lout,
193 : & solout,iout,idid,nmax,uround,meth,nstiff,safe,beta,fac1,fac2,
194 : & work(iek1),work(iek2),work(iek3),work(iek4),work(iek5),
195 : & work(iek6),work(iek7),work(iek8),work(iek9),work(iek10),
196 : & work(iey1),work(ieco),iwork(icomp),nrdens,lrpar,rpar,lipar,ipar,
197 1 : & nfcn,nstep,naccpt,nrejct)
198 1 : work(7)=h
199 1 : iwork(17)=nfcn
200 1 : iwork(18)=nstep
201 1 : iwork(19)=naccpt
202 1 : iwork(20)=nrejct
203 : ! ----------- return -----------
204 1 : return
205 : end subroutine do_dop853
206 : !
207 :
208 : !
209 : ! ----- ... and here is the core integrator ----------
210 : !
211 1 : subroutine dp86co(n,fcn,x,y,xend,hmax,h,rtol,atol,itol,lout,
212 : & solout,iout,idid,nmax,uround,meth,nstiff,safe,beta,fac1,fac2,
213 1 : & k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,y1,rwork,icomp,nrd,lrpar,rpar,lipar,ipar,
214 : & nfcn,nstep,naccpt,nrejct)
215 : ! ----------------------------------------------------------
216 : ! core integrator for dop853
217 : ! parameters same as in dop853 with workspace added
218 : ! ----------------------------------------------------------
219 : ! declarations
220 : ! ----------------------------------------------------------
221 : implicit real(dp) (a-h,o-z)
222 : integer :: n, itol, lout, iout, idid, nmax, meth
223 : parameter (
224 : & c2 = 0.526001519587677318785587544488d-01,
225 : & c3 = 0.789002279381515978178381316732d-01,
226 : & c4 = 0.118350341907227396726757197510d+00,
227 : & c5 = 0.281649658092772603273242802490d+00,
228 : & c6 = 0.333333333333333333333333333333d+00,
229 : & c7 = 0.25d+00,
230 : & c8 = 0.307692307692307692307692307692d+00,
231 : & c9 = 0.651282051282051282051282051282d+00,
232 : & c10 = 0.6d+00,
233 : & c11 = 0.857142857142857142857142857142d+00,
234 : & c14 = 0.1d+00,
235 : & c15 = 0.2d+00,
236 : & c16 = 0.777777777777777777777777777778d+00)
237 : parameter (
238 : & b1 = 5.42937341165687622380535766363d-2,
239 : & b6 = 4.45031289275240888144113950566d0,
240 : & b7 = 1.89151789931450038304281599044d0,
241 : & b8 = -5.8012039600105847814672114227d0,
242 : & b9 = 3.1116436695781989440891606237d-1,
243 : & b10 = -1.52160949662516078556178806805d-1,
244 : & b11 = 2.01365400804030348374776537501d-1,
245 : & b12 = 4.47106157277725905176885569043d-2)
246 : parameter (
247 : & bhh1 = 0.244094488188976377952755905512d+00,
248 : & bhh2 = 0.733846688281611857341361741547d+00,
249 : & bhh3 = 0.220588235294117647058823529412d-01)
250 : parameter (
251 : & er 1 = 0.1312004499419488073250102996d-01,
252 : & er 6 = -0.1225156446376204440720569753d+01,
253 : & er 7 = -0.4957589496572501915214079952d+00,
254 : & er 8 = 0.1664377182454986536961530415d+01,
255 : & er 9 = -0.3503288487499736816886487290d+00,
256 : & er10 = 0.3341791187130174790297318841d+00,
257 : & er11 = 0.8192320648511571246570742613d-01,
258 : & er12 = -0.2235530786388629525884427845d-01)
259 : parameter (
260 : & a21 = 5.26001519587677318785587544488d-2,
261 : & a31 = 1.97250569845378994544595329183d-2,
262 : & a32 = 5.91751709536136983633785987549d-2,
263 : & a41 = 2.95875854768068491816892993775d-2,
264 : & a43 = 8.87627564304205475450678981324d-2,
265 : & a51 = 2.41365134159266685502369798665d-1,
266 : & a53 = -8.84549479328286085344864962717d-1,
267 : & a54 = 9.24834003261792003115737966543d-1,
268 : & a61 = 3.7037037037037037037037037037d-2,
269 : & a64 = 1.70828608729473871279604482173d-1,
270 : & a65 = 1.25467687566822425016691814123d-1,
271 : & a71 = 3.7109375d-2,
272 : & a74 = 1.70252211019544039314978060272d-1,
273 : & a75 = 6.02165389804559606850219397283d-2,
274 : & a76 = -1.7578125d-2)
275 : parameter (
276 : & a81 = 3.70920001185047927108779319836d-2,
277 : & a84 = 1.70383925712239993810214054705d-1,
278 : & a85 = 1.07262030446373284651809199168d-1,
279 : & a86 = -1.53194377486244017527936158236d-2,
280 : & a87 = 8.27378916381402288758473766002d-3,
281 : & a91 = 6.24110958716075717114429577812d-1,
282 : & a94 = -3.36089262944694129406857109825d0,
283 : & a95 = -8.68219346841726006818189891453d-1,
284 : & a96 = 2.75920996994467083049415600797d1,
285 : & a97 = 2.01540675504778934086186788979d1,
286 : & a98 = -4.34898841810699588477366255144d1,
287 : & a101 = 4.77662536438264365890433908527d-1,
288 : & a104 = -2.48811461997166764192642586468d0,
289 : & a105 = -5.90290826836842996371446475743d-1,
290 : & a106 = 2.12300514481811942347288949897d1,
291 : & a107 = 1.52792336328824235832596922938d1,
292 : & a108 = -3.32882109689848629194453265587d1,
293 : & a109 = -2.03312017085086261358222928593d-2)
294 : parameter (
295 : & a111 = -9.3714243008598732571704021658d-1,
296 : & a114 = 5.18637242884406370830023853209d0,
297 : & a115 = 1.09143734899672957818500254654d0,
298 : & a116 = -8.14978701074692612513997267357d0,
299 : & a117 = -1.85200656599969598641566180701d1,
300 : & a118 = 2.27394870993505042818970056734d1,
301 : & a119 = 2.49360555267965238987089396762d0,
302 : & a1110 = -3.0467644718982195003823669022d0,
303 : & a121 = 2.27331014751653820792359768449d0,
304 : & a124 = -1.05344954667372501984066689879d1,
305 : & a125 = -2.00087205822486249909675718444d0,
306 : & a126 = -1.79589318631187989172765950534d1,
307 : & a127 = 2.79488845294199600508499808837d1,
308 : & a128 = -2.85899827713502369474065508674d0,
309 : & a129 = -8.87285693353062954433549289258d0,
310 : & a1210 = 1.23605671757943030647266201528d1,
311 : & a1211 = 6.43392746015763530355970484046d-1)
312 : parameter (
313 : & a141 = 5.61675022830479523392909219681d-2,
314 : & a147 = 2.53500210216624811088794765333d-1,
315 : & a148 = -2.46239037470802489917441475441d-1,
316 : & a149 = -1.24191423263816360469010140626d-1,
317 : & a1410 = 1.5329179827876569731206322685d-1,
318 : & a1411 = 8.20105229563468988491666602057d-3,
319 : & a1412 = 7.56789766054569976138603589584d-3,
320 : & a1413 = -8.298d-3)
321 : parameter (
322 : & a151 = 3.18346481635021405060768473261d-2,
323 : & a156 = 2.83009096723667755288322961402d-2,
324 : & a157 = 5.35419883074385676223797384372d-2,
325 : & a158 = -5.49237485713909884646569340306d-2,
326 : & a1511 = -1.08347328697249322858509316994d-4,
327 : & a1512 = 3.82571090835658412954920192323d-4,
328 : & a1513 = -3.40465008687404560802977114492d-4,
329 : & a1514 = 1.41312443674632500278074618366d-1,
330 : & a161 = -4.28896301583791923408573538692d-1,
331 : & a166 = -4.69762141536116384314449447206d0,
332 : & a167 = 7.68342119606259904184240953878d0,
333 : & a168 = 4.06898981839711007970213554331d0,
334 : & a169 = 3.56727187455281109270669543021d-1,
335 : & a1613 = -1.39902416515901462129418009734d-3,
336 : & a1614 = 2.9475147891527723389556272149d0,
337 : & a1615 = -9.15095847217987001081870187138d0)
338 : parameter (
339 : & d41 = -0.84289382761090128651353491142d+01,
340 : & d46 = 0.56671495351937776962531783590d+00,
341 : & d47 = -0.30689499459498916912797304727d+01,
342 : & d48 = 0.23846676565120698287728149680d+01,
343 : & d49 = 0.21170345824450282767155149946d+01,
344 : & d410 = -0.87139158377797299206789907490d+00,
345 : & d411 = 0.22404374302607882758541771650d+01,
346 : & d412 = 0.63157877876946881815570249290d+00,
347 : & d413 = -0.88990336451333310820698117400d-01,
348 : & d414 = 0.18148505520854727256656404962d+02,
349 : & d415 = -0.91946323924783554000451984436d+01,
350 : & d416 = -0.44360363875948939664310572000d+01)
351 : parameter (
352 : & d51 = 0.10427508642579134603413151009d+02,
353 : & d56 = 0.24228349177525818288430175319d+03,
354 : & d57 = 0.16520045171727028198505394887d+03,
355 : & d58 = -0.37454675472269020279518312152d+03,
356 : & d59 = -0.22113666853125306036270938578d+02,
357 : & d510 = 0.77334326684722638389603898808d+01,
358 : & d511 = -0.30674084731089398182061213626d+02,
359 : & d512 = -0.93321305264302278729567221706d+01,
360 : & d513 = 0.15697238121770843886131091075d+02,
361 : & d514 = -0.31139403219565177677282850411d+02,
362 : & d515 = -0.93529243588444783865713862664d+01,
363 : & d516 = 0.35816841486394083752465898540d+02)
364 : parameter (
365 : & d61 = 0.19985053242002433820987653617d+02,
366 : & d66 = -0.38703730874935176555105901742d+03,
367 : & d67 = -0.18917813819516756882830838328d+03,
368 : & d68 = 0.52780815920542364900561016686d+03,
369 : & d69 = -0.11573902539959630126141871134d+02,
370 : & d610 = 0.68812326946963000169666922661d+01,
371 : & d611 = -0.10006050966910838403183860980d+01,
372 : & d612 = 0.77771377980534432092869265740d+00,
373 : & d613 = -0.27782057523535084065932004339d+01,
374 : & d614 = -0.60196695231264120758267380846d+02,
375 : & d615 = 0.84320405506677161018159903784d+02,
376 : & d616 = 0.11992291136182789328035130030d+02)
377 : parameter (
378 : & d71 = -0.25693933462703749003312586129d+02,
379 : & d76 = -0.15418974869023643374053993627d+03,
380 : & d77 = -0.23152937917604549567536039109d+03,
381 : & d78 = 0.35763911791061412378285349910d+03,
382 : & d79 = 0.93405324183624310003907691704d+02,
383 : & d710 = -0.37458323136451633156875139351d+02,
384 : & d711 = 0.10409964950896230045147246184d+03,
385 : & d712 = 0.29840293426660503123344363579d+02,
386 : & d713 = -0.43533456590011143754432175058d+02,
387 : & d714 = 0.96324553959188282948394950600d+02,
388 : & d715 = -0.39177261675615439165231486172d+02,
389 : & d716 = -0.14972683625798562581422125276d+03)
390 : real(dp) :: y(n),y1(n),k1(n),k2(n),k3(n),k4(n),k5(n),k6(n)
391 : real(dp) :: k7(n),k8(n),k9(n),k10(n),atol(*),rtol(*)
392 2 : dimension icomp(nrd),iwork(nrd+1)
393 : logical :: reject,last
394 : integer :: ierr
395 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
396 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
397 :
398 : interface
399 : #include "num_fcn.dek"
400 : #include "num_solout.dek"
401 : end interface
402 :
403 : real(dp), target :: rwork(3+8*nrd)
404 1 : real(dp), pointer :: cont(:)
405 1 : cont => rwork(3:3+8*nrd)
406 :
407 : ! *** *** *** *** *** *** ***
408 : ! initialisations
409 : ! *** *** *** *** *** *** ***
410 1 : facold=1.d-4
411 1 : expo1=1.d0/8.d0-beta*0.2d0
412 1 : facc1=1.d0/fac1
413 1 : facc2=1.d0/fac2
414 1 : posneg=sign(1.d0,xend-x)
415 : ! --- initial preparations
416 2 : atoli=atol(1)
417 2 : rtoli=rtol(1)
418 1 : last=.false.
419 2 : hlamb=0.d0
420 1 : iasti=0
421 1 : nonsti=0
422 : ierr=0
423 1 : call fcn(n,x,h,y,k1,lrpar,rpar,lipar,ipar,ierr)
424 1 : if (ierr /= 0) then; idid=-5; return; end if
425 1 : hmax=abs(hmax)
426 1 : iord=8
427 1 : if (h == 0.d0) h=hinit(n,fcn,x,y,xend,posneg,k1,k2,k3,iord,
428 0 : & hmax,atol,rtol,itol,lrpar,rpar,lipar,ipar,ierr)
429 1 : if (ierr /= 0) then; idid=-5; return; end if
430 1 : nfcn=nfcn+2
431 1 : reject=.false.
432 2 : xold=x
433 1 : irtrn=1
434 1 : if (iout >= 1) then
435 2 : hout=1.d0
436 1 : rwork(1) = xold
437 1 : rwork(2) = hout
438 1 : iwork(1) = nrd
439 3 : iwork(2:nrd+1) = icomp(1:nrd)
440 :
441 3 : do j=1,nrd
442 2 : i=icomp(j)
443 2 : cont(j)=y(i)
444 2 : cont(j+nrd)=0
445 2 : cont(j+nrd*2)=0
446 2 : cont(j+nrd*3)=0
447 2 : cont(j+nrd*4)=0
448 2 : cont(j+nrd*5)=0
449 2 : cont(j+nrd*6)=0
450 3 : cont(j+nrd*7)=0
451 : end do
452 :
453 1 : call solout(naccpt+1,xold,x,n,y,rwork,iwork,contd8,lrpar,rpar,lipar,ipar,irtrn)
454 1 : if (irtrn <0) GOTO 79
455 : end if
456 : ! --- basic integration step
457 : 1 continue
458 663 : if (nstep >nmax) GOTO 78
459 663 : if (0.1d0*abs(h) <= abs(x)*uround)GOTO 77
460 663 : if ((x+1.01d0*h-xend)*posneg >0.d0) then
461 1 : h=xend-x
462 1 : last=.true.
463 : end if
464 663 : nstep=nstep+1
465 : ! --- the twelve stages
466 663 : if (irtrn >= 2) then
467 0 : call fcn(n,x,h,y,k1,lrpar,rpar,lipar,ipar,ierr)
468 0 : if (ierr /= 0) then; hnew=h/facc1; h=hnew; GOTO 1; end if
469 : end if
470 1989 : do i=1,n
471 1989 : y1(i)=y(i)+h*a21*k1(i)
472 : end do
473 663 : call fcn(n,x+c2*h,h,y1,k2,lrpar,rpar,lipar,ipar,ierr)
474 663 : if (ierr /= 0) then; hnew=h/facc1; h=hnew; GOTO 1; end if
475 1989 : do i=1,n
476 1989 : y1(i)=y(i)+h*(a31*k1(i)+a32*k2(i))
477 : end do
478 663 : call fcn(n,x+c3*h,h,y1,k3,lrpar,rpar,lipar,ipar,ierr)
479 663 : if (ierr /= 0) then; hnew=h/facc1; h=hnew; GOTO 1; end if
480 1989 : do i=1,n
481 1989 : y1(i)=y(i)+h*(a41*k1(i)+a43*k3(i))
482 : end do
483 663 : call fcn(n,x+c4*h,h,y1,k4,lrpar,rpar,lipar,ipar,ierr)
484 663 : if (ierr /= 0) then; hnew=h/facc1; h=hnew; GOTO 1; end if
485 1989 : do i=1,n
486 1989 : y1(i)=y(i)+h*(a51*k1(i)+a53*k3(i)+a54*k4(i))
487 : end do
488 663 : call fcn(n,x+c5*h,h,y1,k5,lrpar,rpar,lipar,ipar,ierr)
489 663 : if (ierr /= 0) then; hnew=h/facc1; h=hnew; GOTO 1; end if
490 1989 : do i=1,n
491 1989 : y1(i)=y(i)+h*(a61*k1(i)+a64*k4(i)+a65*k5(i))
492 : end do
493 663 : call fcn(n,x+c6*h,h,y1,k6,lrpar,rpar,lipar,ipar,ierr)
494 663 : if (ierr /= 0) then; hnew=h/facc1; h=hnew; GOTO 1; end if
495 1989 : do i=1,n
496 1989 : y1(i)=y(i)+h*(a71*k1(i)+a74*k4(i)+a75*k5(i)+a76*k6(i))
497 : end do
498 663 : call fcn(n,x+c7*h,h,y1,k7,lrpar,rpar,lipar,ipar,ierr)
499 663 : if (ierr /= 0) then; hnew=h/facc1; h=hnew; GOTO 1; end if
500 1989 : do i=1,n
501 1989 : y1(i)=y(i)+h*(a81*k1(i)+a84*k4(i)+a85*k5(i)+a86*k6(i)+a87*k7(i))
502 : end do
503 663 : call fcn(n,x+c8*h,h,y1,k8,lrpar,rpar,lipar,ipar,ierr)
504 663 : if (ierr /= 0) then; hnew=h/facc1; h=hnew; GOTO 1; end if
505 1989 : do i=1,n
506 : y1(i)=y(i)+h*(a91*k1(i)+a94*k4(i)+a95*k5(i)+a96*k6(i)+a97*k7(i)
507 1989 : & +a98*k8(i))
508 : end do
509 663 : call fcn(n,x+c9*h,h,y1,k9,lrpar,rpar,lipar,ipar,ierr)
510 663 : if (ierr /= 0) then; hnew=h/facc1; h=hnew; GOTO 1; end if
511 1989 : do i=1,n
512 : y1(i)=y(i)+h*(a101*k1(i)+a104*k4(i)+a105*k5(i)+a106*k6(i)
513 1989 : & +a107*k7(i)+a108*k8(i)+a109*k9(i))
514 : end do
515 663 : call fcn(n,x+c10*h,h,y1,k10,lrpar,rpar,lipar,ipar,ierr)
516 663 : if (ierr /= 0) then; hnew=h/facc1; h=hnew; GOTO 1; end if
517 1989 : do i=1,n
518 : y1(i)=y(i)+h*(a111*k1(i)+a114*k4(i)+a115*k5(i)+a116*k6(i)
519 1989 : & +a117*k7(i)+a118*k8(i)+a119*k9(i)+a1110*k10(i))
520 : end do
521 663 : call fcn(n,x+c11*h,h,y1,k2,lrpar,rpar,lipar,ipar,ierr)
522 663 : if (ierr /= 0) then; hnew=h/facc1; h=hnew; GOTO 1; end if
523 664 : xph=x+h
524 1989 : do i=1,n
525 : y1(i)=y(i)+h*(a121*k1(i)+a124*k4(i)+a125*k5(i)+a126*k6(i)
526 1989 : & +a127*k7(i)+a128*k8(i)+a129*k9(i)+a1210*k10(i)+a1211*k2(i))
527 : end do
528 663 : call fcn(n,xph,h,y1,k3,lrpar,rpar,lipar,ipar,ierr)
529 663 : if (ierr /= 0) then; hnew=h/facc1; h=hnew; GOTO 1; end if
530 663 : nfcn=nfcn+11
531 1989 : do i=1,n
532 : k4(i)=b1*k1(i)+b6*k6(i)+b7*k7(i)+b8*k8(i)+b9*k9(i)
533 1326 : & +b10*k10(i)+b11*k2(i)+b12*k3(i)
534 1989 : k5(i)=y(i)+h*k4(i)
535 : end do
536 : ! --- error estimation
537 664 : err=0.d0
538 664 : err2=0.d0
539 663 : if (itol == 0) then
540 1989 : do i=1,n
541 1327 : sk=atoli+rtoli*max(abs(y(i)),abs(k5(i)))
542 1327 : erri=k4(i)-bhh1*k1(i)-bhh2*k9(i)-bhh3*k3(i)
543 1326 : err2=err2+(erri/sk)**2
544 : erri=er1*k1(i)+er6*k6(i)+er7*k7(i)+er8*k8(i)+er9*k9(i)
545 1326 : & +er10*k10(i)+er11*k2(i)+er12*k3(i)
546 1989 : err=err+(erri/sk)**2
547 : end do
548 : else
549 0 : do i=1,n
550 0 : sk=atol(i)+rtol(i)*max(abs(y(i)),abs(k5(i)))
551 0 : erri=k4(i)-bhh1*k1(i)-bhh2*k9(i)-bhh3*k3(i)
552 0 : err2=err2+(erri/sk)**2
553 : erri=er1*k1(i)+er6*k6(i)+er7*k7(i)+er8*k8(i)+er9*k9(i)
554 0 : & +er10*k10(i)+er11*k2(i)+er12*k3(i)
555 0 : err=err+(erri/sk)**2
556 : end do
557 : end if
558 664 : deno=err+0.01d0*err2
559 663 : if (deno <= 0.d0) deno=1.d0
560 663 : err=abs(h)*err*sqrt(1.d0/(n*deno))
561 : ! --- computation of hnew
562 664 : fac11=pow(err,expo1)
563 : ! --- lund-stabilization
564 664 : fac=fac11/pow(facold,beta)
565 : ! --- we require fac1 <= hnew/h <= fac2
566 663 : fac=max(facc2,min(facc1,fac/safe))
567 663 : hnew=h/fac
568 663 : if(err <= 1.d0)then
569 : ! --- step is accepted
570 627 : facold=max(err,1.0d-4)
571 627 : naccpt=naccpt+1
572 627 : call fcn(n,xph,h,k5,k4,lrpar,rpar,lipar,ipar,ierr)
573 627 : if (ierr /= 0) then; hnew=h/facc1; h=hnew; GOTO 1; end if
574 627 : nfcn=nfcn+1
575 : ! ------- stiffness detection
576 627 : if (mod(naccpt,nstiff) == 0.or.iasti >0) then
577 628 : stnum=0.d0
578 628 : stden=0.d0
579 1881 : do i=1,n
580 1254 : stnum=stnum+(k4(i)-k3(i))**2
581 1881 : stden=stden+(k5(i)-y1(i))**2
582 : end do
583 627 : if (stden >0.d0) hlamb=abs(h)*sqrt(stnum/stden)
584 627 : if (hlamb >6.1d0) then
585 550 : nonsti=0
586 550 : iasti=iasti+1
587 550 : if (iasti == 15) then
588 3 : if (lout >0) write (lout,*)
589 0 : & ' the problem seems to become stiff at x = ',x
590 3 : if (lout <0) GOTO 76
591 : end if
592 : else
593 77 : nonsti=nonsti+1
594 77 : if (nonsti == 6) iasti=0
595 : end if
596 : end if
597 : ! ------- final preparation for dense output
598 627 : if (iout >= 2) then
599 : ! ---- save the first function evaluations
600 1881 : do j=1,nrd
601 1254 : i=icomp(j)
602 1254 : cont(j)=y(i)
603 1255 : ydiff=k5(i)-y(i)
604 1254 : cont(j+nrd)=ydiff
605 1255 : bspl=h*k1(i)-ydiff
606 1254 : cont(j+nrd*2)=bspl
607 1254 : cont(j+nrd*3)=ydiff-h*k4(i)-bspl
608 : cont(j+nrd*4)=d41*k1(i)+d46*k6(i)+d47*k7(i)+d48*k8(i)
609 1254 : & +d49*k9(i)+d410*k10(i)+d411*k2(i)+d412*k3(i)
610 : cont(j+nrd*5)=d51*k1(i)+d56*k6(i)+d57*k7(i)+d58*k8(i)
611 1254 : & +d59*k9(i)+d510*k10(i)+d511*k2(i)+d512*k3(i)
612 : cont(j+nrd*6)=d61*k1(i)+d66*k6(i)+d67*k7(i)+d68*k8(i)
613 1254 : & +d69*k9(i)+d610*k10(i)+d611*k2(i)+d612*k3(i)
614 : cont(j+nrd*7)=d71*k1(i)+d76*k6(i)+d77*k7(i)+d78*k8(i)
615 1881 : & +d79*k9(i)+d710*k10(i)+d711*k2(i)+d712*k3(i)
616 : end do
617 : ! --- the next three function evaluations
618 1881 : do i=1,n
619 : y1(i)=y(i)+h*(a141*k1(i)+a147*k7(i)+a148*k8(i)
620 : & +a149*k9(i)+a1410*k10(i)+a1411*k2(i)+a1412*k3(i)
621 1881 : & +a1413*k4(i))
622 : end do
623 627 : call fcn(n,x+c14*h,h,y1,k10,lrpar,rpar,lipar,ipar,ierr)
624 627 : if (ierr /= 0) then; hnew=h/facc1; h=hnew; GOTO 1; end if
625 1881 : do i=1,n
626 : y1(i)=y(i)+h*(a151*k1(i)+a156*k6(i)+a157*k7(i)
627 : & +a158*k8(i)+a1511*k2(i)+a1512*k3(i)+a1513*k4(i)
628 1881 : & +a1514*k10(i))
629 : end do
630 627 : call fcn(n,x+c15*h,h,y1,k2,lrpar,rpar,lipar,ipar,ierr)
631 627 : if (ierr /= 0) then; hnew=h/facc1; h=hnew; GOTO 1; end if
632 1881 : do i=1,n
633 : y1(i)=y(i)+h*(a161*k1(i)+a166*k6(i)+a167*k7(i)
634 : & +a168*k8(i)+a169*k9(i)+a1613*k4(i)+a1614*k10(i)
635 1881 : & +a1615*k2(i))
636 : end do
637 627 : call fcn(n,x+c16*h,h,y1,k3,lrpar,rpar,lipar,ipar,ierr)
638 627 : if (ierr /= 0) then; hnew=h/facc1; h=hnew; GOTO 1; end if
639 627 : nfcn=nfcn+3
640 : ! --- final preparation
641 1881 : do j=1,nrd
642 1254 : i=icomp(j)
643 : cont(j+nrd*4)=h*(cont(j+nrd*4)+d413*k4(i)+d414*k10(i)
644 1254 : & +d415*k2(i)+d416*k3(i))
645 : cont(j+nrd*5)=h*(cont(j+nrd*5)+d513*k4(i)+d514*k10(i)
646 1254 : & +d515*k2(i)+d516*k3(i))
647 : cont(j+nrd*6)=h*(cont(j+nrd*6)+d613*k4(i)+d614*k10(i)
648 1254 : & +d615*k2(i)+d616*k3(i))
649 : cont(j+nrd*7)=h*(cont(j+nrd*7)+d713*k4(i)+d714*k10(i)
650 1881 : & +d715*k2(i)+d716*k3(i))
651 : end do
652 627 : hout=h
653 : end if
654 1881 : do i=1,n
655 1254 : k1(i)=k4(i)
656 1881 : y(i)=k5(i)
657 : end do
658 627 : xold=x
659 627 : x=xph
660 627 : if (iout >= 1) then
661 627 : rwork(1) = xold
662 627 : rwork(2) = hout
663 627 : iwork(1) = nrd
664 1881 : iwork(2:nrd+1) = icomp(1:nrd)
665 627 : call solout(naccpt+1,xold,x,n,y,rwork,iwork,contd8,lrpar,rpar,lipar,ipar,irtrn)
666 627 : if (irtrn <0) GOTO 79
667 : end if
668 : ! ------- normal exit
669 627 : if (last) then
670 1 : h=hnew
671 1 : idid=1
672 1 : return
673 : end if
674 626 : if(abs(hnew) >hmax)hnew=posneg*hmax
675 626 : if(reject)hnew=posneg*min(abs(hnew),abs(h))
676 : reject=.false.
677 : else
678 : ! --- step is rejected
679 36 : hnew=h/min(facc1,fac11/safe)
680 36 : reject=.true.
681 36 : if(naccpt >= 1)nrejct=nrejct+1
682 : last=.false.
683 : end if
684 662 : h=hnew
685 662 : GOTO 1
686 : ! --- fail exit
687 : 76 continue
688 0 : idid=-4
689 0 : return
690 : 77 continue
691 0 : if (lout >0) write(lout,979) x
692 0 : if (lout >0) write(lout,*)' step size too small, h=',h
693 0 : idid=-3
694 0 : return
695 : 78 continue
696 0 : if (lout >0) write(lout,979) x
697 0 : if (lout >0) write(lout,*)
698 0 : & ' more than nmax =',nmax,'steps are needed'
699 0 : idid=-2
700 0 : return
701 : 79 continue
702 : !if (lout >0) write(lout,979) x
703 : 979 format(' exit of dop853 at x=',e18.4)
704 0 : idid=2
705 0 : return
706 1 : end subroutine dp86co
707 : !
708 0 : function hinit(n,fcn,x,y,xend,posneg,f0,f1,y1,iord,
709 : & hmax,atol,rtol,itol,lrpar,rpar,lipar,ipar,ierr)
710 : ! ----------------------------------------------------------
711 : ! ---- computation of an initial step size guess
712 : ! ----------------------------------------------------------
713 : implicit real(dp) (a-h,o-z)
714 :
715 : real(dp), intent(inout) :: x
716 : real(dp), intent(inout), dimension(:) :: y, f0, f1 ! (n)
717 : real(dp), intent(in) :: xend
718 : real(dp), intent(in) :: rtol(*)
719 : real(dp), intent(in) :: atol(*)
720 : integer, intent(in) :: itol
721 : integer, intent(in) :: lrpar, lipar
722 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
723 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
724 : integer, intent(out) :: ierr
725 :
726 : dimension y1(n)
727 :
728 : interface
729 : #include "num_fcn.dek"
730 : end interface
731 :
732 : ! ---- compute a first guess for explicit euler as
733 : ! ---- h = 0.01 * norm (y0) / norm (f0)
734 : ! ---- the increment for explicit euler is small
735 : ! ---- compared to the solution
736 0 : dnf=0.0d0
737 0 : dny=0.0d0
738 0 : atoli=atol(1)
739 0 : rtoli=rtol(1)
740 0 : if (itol == 0) then
741 0 : do i=1,n
742 0 : sk=atoli+rtoli*abs(y(i))
743 0 : dnf=dnf+(f0(i)/sk)**2
744 0 : dny=dny+(y(i)/sk)**2
745 : end do
746 : else
747 0 : do i=1,n
748 0 : sk=atol(i)+rtol(i)*abs(y(i))
749 0 : dnf=dnf+(f0(i)/sk)**2
750 0 : dny=dny+(y(i)/sk)**2
751 : end do
752 : end if
753 0 : if (dnf <= 1.d-10.or.dny <= 1.d-10) then
754 0 : h=1.0d-6
755 : else
756 0 : h=sqrt(dny/dnf)*0.01d0
757 : end if
758 0 : h=min(h,hmax)
759 0 : h=sign(h,posneg)
760 : ! ---- perform an explicit euler step
761 0 : do i=1,n
762 0 : y1(i)=y(i)+h*f0(i)
763 : end do
764 0 : call fcn(n,x+h,h,y1,f1,lrpar,rpar,lipar,ipar,ierr)
765 0 : if (ierr /= 0) then; idid=-5; return; end if
766 : ! ---- estimate the second derivative of the solution
767 0 : der2=0.0d0
768 0 : if (itol == 0) then
769 0 : do i=1,n
770 0 : sk=atoli+rtoli*abs(y(i))
771 0 : der2=der2+((f1(i)-f0(i))/sk)**2
772 : end do
773 : else
774 0 : do i=1,n
775 0 : sk=atol(i)+rtol(i)*abs(y(i))
776 0 : der2=der2+((f1(i)-f0(i))/sk)**2
777 : end do
778 : end if
779 0 : der2=sqrt(der2)/h
780 : ! ---- step size is computed such that
781 : ! ---- h**iord * max ( norm (f0), norm (der2)) = 0.01
782 0 : der12=max(abs(der2),sqrt(dnf))
783 0 : if (der12 <= 1.d-15) then
784 0 : h1=max(1.0d-6,abs(h)*1.0d-3)
785 : else
786 0 : h1=pow(0.01d0/der12,1.d0/iord)
787 : end if
788 0 : h=min(100*abs(h),h1,hmax)
789 0 : hinit=sign(h,posneg)
790 : return
791 : end function hinit
792 :
793 0 : real(dp) function contd8(ii,x,rwork,iwork,ierr)
794 : ! ----------------------------------------------------------
795 : ! this function can be used for continuous output in connection
796 : ! with the output-subroutine for dop853. it provides an
797 : ! approximation to the ii-th component of the solution at x.
798 : ! ----------------------------------------------------------
799 : integer, intent(in) :: ii ! result is interpolated approximation of y(i) at x=s.
800 : real(dp), intent(in) :: x ! interpolation x value (between xold and x).
801 : real(dp), intent(inout), target :: rwork(*)
802 : integer, intent(inout), target :: iwork(*)
803 : integer, intent(out) :: ierr
804 :
805 0 : real(dp) :: xold, h, s, s1
806 : integer :: nd, i, j
807 0 : real(dp), pointer :: con(:)
808 0 : integer, pointer :: icomp(:)
809 :
810 0 : nd = iwork(1)
811 0 : icomp => iwork(2:nd+1)
812 0 : xold = rwork(1)
813 0 : h = rwork(2)
814 0 : con => rwork(3:3+8*nd)
815 :
816 : ! ----- compute place of ii-th component
817 0 : i=0
818 0 : do j=1,nd
819 0 : if (icomp(j) == ii) then
820 : i=j; exit
821 : end if
822 : end do
823 0 : if (i == 0) then
824 0 : contd8 = 0
825 0 : ierr = -1
826 0 : return
827 : end if
828 0 : ierr=0
829 0 : s=(x-xold)/h
830 0 : s1=1.d0-s
831 0 : conpar=con(i+nd*4)+s*(con(i+nd*5)+s1*(con(i+nd*6)+s*con(i+nd*7)))
832 0 : contd8=con(i)+s*(con(i+nd)+s1*(con(i+nd*2)+s*(con(i+nd*3)+s1*conpar)))
833 :
834 0 : end function contd8
835 :
836 : end module mod_dop853
|
