MPI-AMRVAC 3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
Loading...
Searching...
No Matches
mod_odeint.t
Go to the documentation of this file.
1!> This module packages odeint from numerical recipes.
3 use mod_comm_lib, only: mpistop
4 implicit none
5 double precision :: tiny
6 parameter(tiny=1.d-30)
7 integer :: maxstp, nmax, kmaxx
8 parameter(maxstp=10000,nmax=50,kmaxx=200)
9
10contains
11
12 subroutine odeint(ystart,nvar,x1,x2,eps,h1,hmin,nok,nbad,derivs,rkqs,ierror)
13
14 INTEGER nbad,nok,nvar
15 double precision eps,h1,hmin,x1,x2,ystart(nvar)
16 double precision h,hdid,hnext,x,xsav,dydx(NMAX),y(NMAX),yscal(NMAX)
17 double precision dxsav, yp(NMAX,KMAXX), xp(KMAXX)
18 integer kmax, kount
19 INTEGER i,nstp
20 ! Can be 0: success, 1: hmin too small, 2: MAXSTP exceeded
21 integer, intent(out) :: ierror
22
23 EXTERNAL derivs,rkqs
24
25 COMMON /path/ kmax,kount,dxsav,xp,yp
26
27 x = x1
28 h = sign(h1,x2-x1)
29 nok = 0
30 nbad = 0
31 kount = 0
32 xsav = 0.d0
33
34 do i=1,nvar
35 y(i)=ystart(i)
36 end do
37
38 if (kmax.gt.0) xsav=x-2.d0*dxsav
39
40 do nstp=1,maxstp
41 call derivs(x,y,dydx)
42 do i=1,nvar
43 yscal(i)=abs(y(i))+abs(h*dydx(i))+tiny
44 end do
45
46 if(kmax.gt.0)then
47 if(abs(x-xsav).gt.abs(dxsav)) then
48 if(kount.lt.kmax-1)then
49 kount=kount+1
50 xp(kount)=x
51 do i=1,nvar
52 yp(i,kount)=y(i)
53 end do
54 xsav=x
55 end if
56 end if
57 end if
58
59 if((x+h-x2)*(x+h-x1).gt.0.d0) h=x2-x
60
61 if (any(y(1:nvar) .ne. y(1:nvar)) .or. any(dydx(1:nvar) .ne. dydx(1:nvar)) .or. x .ne. x) then
62 write(*,*) "ODEINT BEFORE CALL RKQS: NaN in x, dydx, or y!"
63 write(*,*) "x",x
64 write(*,*) "y",y
65 write(*,*) "dydx",dydx
66 write(*,*) "exiting..."
67 ierror = 2
68 return
69 end if
70
71 call rkqs(y,dydx,nvar,x,h,eps,yscal,hdid,hnext,derivs)
72
73 if (any(y(1:nvar) .ne. y(1:nvar)) .or. any(dydx(1:nvar) .ne. dydx(1:nvar)) .or. x .ne. x) then
74 write(*,*) "ODEINT AFTER CALL RKQS: NaN in x, dydx, or y!"
75 write(*,*) "x",x
76 write(*,*) "y",y
77 write(*,*) "dydx",dydx
78 write(*,*) "exiting..."
79 ierror = 2
80 return
81 end if
82
83 if(hdid.eq.h)then
84 nok=nok+1
85 else
86 nbad=nbad+1
87 end if
88
89 if((x-x2)*(x2-x1).ge.0.d0)then
90 do i=1,nvar
91 ystart(i)=y(i)
92 end do
93 if(kmax.ne.0)then
94 kount=kount+1
95 xp(kount)=x
96 do i=1,nvar
97 yp(i,kount)=y(i)
98 end do
99 end if
100
101 ierror = 0
102 return
103 end if
104
105 if(abs(hnext).lt.hmin) then
106 ierror = 1
107 end if
108
109 h=hnext
110 end do
111
112 ierror = 2
113 end subroutine odeint
114
115 SUBROUTINE rkqs(y,dydx,n,x,htry,eps,yscal,hdid,hnext,derivs)
116
117 INTEGER n,NMAX
118 double precision eps,hdid,hnext,htry,x,dydx(n),y(n),yscal(n)
119 EXTERNAL derivs
120 parameter(nmax=50)
121! CU USES derivs,rkck
122 INTEGER i
123 double precision errmax,h,xnew,yerr(NMAX),ytemp(NMAX),SAFETY,PGROW,PSHRNK,ERRCON
124 parameter(safety=0.9d0,pgrow=-.2d0,pshrnk=-.25d0,errcon=1.89d-4)
125
126 h=htry
127
128 if (any(y(1:n) .ne. y(1:n)) .or. any(dydx(1:n) .ne. dydx(1:n)) .or. x .ne. x) then
129 write(*,*) "RKQS BEFORE CALL RKCK: NaN in x, dydx, or y!"
130 write(*,*) "x",x
131 write(*,*) "y",y
132 write(*,*) "dydx",dydx
133 write(*,*) "exiting..."
134 return
135 end if
136
137 1 call rkck(y,dydx,n,x,h,ytemp,yerr,derivs)
138 if (any(y(1:n) .ne. y(1:n)) .or. any(dydx(1:n) .ne. dydx(1:n)) .or. x .ne. x) then
139 write(*,*) "RKQS AFTER CALL RKCK: NaN in x, dydx, or y!"
140 write(*,*) "x",x
141 write(*,*) "y",y
142 write(*,*) "dydx",dydx
143 write(*,*) "exiting..."
144 return
145 end if
146
147 errmax=0.d0
148 do 11 i=1,n
149 errmax=max(errmax,abs(yerr(i)/yscal(i)))
150 11 continue
151 errmax=errmax/eps
152 if(errmax.gt.1.d0)then
153 h=safety*h*(errmax**pshrnk)
154 if(h.lt.0.1d0*h)then
155 h=.1d0*h
156 end if
157 xnew=x+h
158 if(xnew.eq.x) call stop('stepsize underflow in rkqs')
159 goto 1
160 else
161 if(errmax.gt.errcon)then
162 hnext=safety*h*(errmax**pgrow)
163 else
164 hnext=5.d0*h
165 end if
166 hdid=h
167 x=x+h
168 do 12 i=1,n
169 y(i)=ytemp(i)
170 12 continue
171 return
172 end if
173 END subroutine rkqs
174!C (C) Copr. 1986-92 Numerical Recipes Software Vs1&v%1jw#<?4210(9p#.
175
176 SUBROUTINE rkck(y,dydx,n,x,h,yout,yerr,derivs)
177
178 INTEGER n,NMAX
179 double precision h,x,dydx(n),y(n),yerr(n),yout(n)
180 EXTERNAL derivs
181 parameter(nmax=50)
182!CU USES derivs
183 INTEGER i
184 double precision ak2(NMAX),ak3(NMAX),ak4(NMAX),ak5(NMAX),ak6(NMAX),ytemp(NMAX),A2,A3,A4,A5,A6,B21,B31,B32,B41,B42,B43,B51, &
185 B52,B53,&
186 B54,B61,B62,B63,B64,B65,C1,C3,C4,C6,DC1,DC3,DC4,DC5,DC6
187 parameter(a2=.2d0,a3=.3d0,a4=.6d0,a5=1.d0, &
188 a6=.875d0,b21=.2d0,b31=3.d0/40.d0,&
189 b32=9.d0/40.d0,b41=.3d0,b42=-.9d0,b43=1.2d0,&
190 b51=-11.d0/54.d0,b52=2.5d0,&
191 b53=-70.d0/27.d0,b54=35.d0/27.d0,b61=1631.d0/55296.d0, &
192 b62=175.d0/512.d0, &
193 b63=575.d0/13824.d0,b64=44275.d0/110592.d0, &
194 b65=253.d0/4096.d0,c1=37.d0/378.d0,&
195 c3=250.d0/621.d0,c4=125.d0/594.d0,c6=512.d0/1771.d0,&
196 dc1=c1-2825.d0/27648.d0,&
197 dc3=c3-18575.d0/48384.d0,dc4=c4-13525.d0/55296.d0,&
198 dc5=-277.d0/14336.d0,&
199 dc6=c6-.25d0)
200 if (any(y(1:n) .ne. y(1:n)) .or. any(dydx(1:n) .ne. dydx(1:n))) then
201 write(*,*) "NaNs IN RKCK, STEP 0!"
202 write(*,*) "y0",y(1:n)
203 write(*,*) "derivs",dydx(1:n)
204 write(*,*) "ABORTING..."
205 call mpistop("")
206 end if
207 do 11 i=1,n
208 ytemp(i)=y(i)+b21*h*dydx(i)
209 11 continue
210 call derivs(x+a2*h,ytemp,ak2)
211 if (any(ytemp(1:n) .ne. ytemp(1:n)) .or. any(ak2(1:n) .ne. ak2(1:n))) then
212 write(*,*) "NaNs IN RKCK, STEP 1!"
213 write(*,*) "y1",ytemp(1:n)
214 write(*,*) "derivs",ak2(1:n)
215 write(*,*) "ABORTING..."
216 call mpistop("")
217 end if
218 do 12 i=1,n
219 ytemp(i)=y(i)+h*(b31*dydx(i)+b32*ak2(i))
220 12 continue
221 call derivs(x+a3*h,ytemp,ak3)
222 if (any(ytemp(1:n) .ne. ytemp(1:n)) .or. any(ak3(1:n) .ne. ak3(1:n))) then
223 write(*,*) "NaNs IN RKCK, STEP 2!"
224 write(*,*) "y2",ytemp(1:n)
225 write(*,*) "derivs",ak3(1:n)
226 write(*,*) "ABORTING..."
227 call mpistop("")
228 end if
229 do 13 i=1,n
230 ytemp(i)=y(i)+h*(b41*dydx(i)+b42*ak2(i)+b43*ak3(i))
231 13 continue
232 call derivs(x+a4*h,ytemp,ak4)
233 if (any(ytemp(1:n) .ne. ytemp(1:n)) .or. any(ak4(1:n) .ne. ak4(1:n))) then
234 write(*,*) "NaNs IN RKCK, STEP 3!"
235 write(*,*) "y3",ytemp(1:n)
236 write(*,*) "derivs",ak4(1:n)
237 write(*,*) "ABORTING..."
238 call mpistop("")
239 end if
240 do 14 i=1,n
241 ytemp(i)=y(i)+h*(b51*dydx(i)+b52*ak2(i)+b53*ak3(i)+b54*ak4(i))
242 14 continue
243 call derivs(x+a5*h,ytemp,ak5)
244 if (any(ytemp(1:n) .ne. ytemp(1:n)) .or. any(ak5(1:n) .ne. ak5(1:n))) then
245 write(*,*) "NaNs IN RKCK, STEP 4!"
246 write(*,*) "y4",ytemp(1:n)
247 write(*,*) "derivs",ak5(1:n)
248 write(*,*) "ABORTING..."
249 call mpistop("")
250 end if
251 do 15 i=1,n
252 ytemp(i)=y(i)+h*(b61*dydx(i)+b62*ak2(i)+b63*ak3(i)+b64*ak4(i)+b65*ak5(i))
253 15 continue
254 call derivs(x+a6*h,ytemp,ak6)
255 if (any(ytemp(1:n) .ne. ytemp(1:n)) .or. any(ak6(1:n) .ne. ak6(1:n))) then
256 write(*,*) "NaNs IN RKCK, STEP 1!"
257 write(*,*) "y15",ytemp(1:n)
258 write(*,*) "derivs",ak6(1:n)
259 write(*,*) "ABORTING..."
260 call mpistop("")
261 end if
262 do 16 i=1,n
263 yout(i)=y(i)+h*(c1*dydx(i)+c3*ak3(i)+c4*ak4(i)+c6*ak6(i))
264 16 continue
265 do 17 i=1,n
266 yerr(i)=h*(dc1*dydx(i)+dc3*ak3(i)+dc4*ak4(i)+dc5*ak5(i)+dc6*ak6(i))
267 17 continue
268 return
269 END subroutine rkck
270!C (C) Copr. 1986-92 Numerical Recipes Software Vs1&v%1jw#<?4210(9p#.
271 subroutine stop(text)
272 character(len=*), intent(in) :: text
273 print*, text
274 stop
275 end subroutine stop
276
277end module mod_odeint
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
This module packages odeint from numerical recipes.
Definition mod_odeint.t:2
subroutine odeint(ystart, nvar, x1, x2, eps, h1, hmin, nok, nbad, derivs, rkqs, ierror)
Definition mod_odeint.t:13
integer nmax
Definition mod_odeint.t:7
subroutine stop(text)
Definition mod_odeint.t:272
integer maxstp
Definition mod_odeint.t:7
integer kmaxx
Definition mod_odeint.t:7
subroutine rkqs(y, dydx, n, x, htry, eps, yscal, hdid, hnext, derivs)
Definition mod_odeint.t:116
subroutine rkck(y, dydx, n, x, h, yout, yerr, derivs)
Definition mod_odeint.t:177
double precision tiny
Definition mod_odeint.t:5