MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
mod_odeint.t
Go to the documentation of this file.
1 !> This module packages odeint from numerical recipes.
2 module mod_odeint
3  use mod_comm_lib, only: mpistop
4  implicit none
5  integer :: maxstp, nmax, kmaxx
6  parameter(maxstp=10000,nmax=50,kmaxx=200)
7  double precision :: tiny
8  parameter(tiny=1.d-30)
9 
10 contains
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  INTEGER i,nstp
17  double precision h,hdid,hnext,x,xsav,dydx(NMAX),y(NMAX),yscal(NMAX)
18  integer kmax, kount
19  double precision dxsav, yp(NMAX,KMAXX), xp(KMAXX)
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 
277 end module mod_odeint
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: mod_comm_lib.t:208
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:5
subroutine stop(text)
Definition: mod_odeint.t:272
integer maxstp
Definition: mod_odeint.t:5
integer kmaxx
Definition: mod_odeint.t:5
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:7