MPI-AMRVAC 3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
Loading...
Searching...
No Matches
mod_tdfluxrope.t
Go to the documentation of this file.
1module mod_tdfluxrope
2 implicit none
3 double precision :: d_td99,l_td99,r_td99,a_td99,q_td99,izero_td99,li_td99
4 double precision :: p_bt_ratio
5
6contains
7
8{^ifthreed
9 subroutine td99(ixI^L,ixO^L,x,B)
10 use mod_global_parameters
11
12 integer, intent(in) :: ixI^L, ixO^L
13 double precision, intent(in) :: x(ixI^S,1:ndim)
14 double precision, intent(inout) :: B(ixI^S,1:ndir)
15
16 b=0.d0
17 call td99_bi (ixi^l,ixo^l,x,b)
18 call td99_bth(ixi^l,ixo^l,x,b)
19 call td99_bq (ixi^l,ixo^l,x,b)
20
21 end subroutine td99
22
23 subroutine td99_bq(ixI^L,ixO^L,x,B)
24 use mod_global_parameters
25
26 integer, intent(in) :: ixI^L, ixO^L
27 double precision, intent(in) :: x(ixI^S,1:ndim)
28 double precision, intent(inout) :: B(ixI^S,1:ndir)
29
30 double precision :: rplus(ixI^S,1:ndir),rminu(ixI^S,1:ndir)
31 double precision :: rplusvalue(ixI^S),rminuvalue(ixI^S)
32 integer :: idir
33
34 rplus(ixo^s,1)=x(ixo^s,1)-l_td99
35 rminu(ixo^s,1)=x(ixo^s,1)+l_td99
36 rplus(ixo^s,2)=x(ixo^s,2)
37 rminu(ixo^s,2)=x(ixo^s,2)
38 rplus(ixo^s,3)=x(ixo^s,3)+d_td99
39 rminu(ixo^s,3)=x(ixo^s,3)+d_td99
40 rplusvalue(ixo^s)=dsqrt(^d&rplus(ixo^s,^d)**2+)
41 rminuvalue(ixo^s)=dsqrt(^d&rminu(ixo^s,^d)**2+)
42 do idir=1,ndir
43 b(ixo^s,idir)=b(ixo^s,idir)+q_td99 &
44 /sqrt(4.0d0*dpi)*(rplus(ixo^s,idir) &
45 /rplusvalue(ixo^s)**3-rminu(ixo^s,idir)/rminuvalue(ixo^s)**3)
46 end do
47
48 end subroutine td99_bq
49
50 subroutine td99_bth(ixI^L,ixO^L,x,B)
51 use mod_global_parameters
52 integer, intent(in) :: ixI^L, ixO^L
53 double precision, intent(in) :: x(ixI^S,1:ndim)
54 double precision, intent(inout) :: B(ixI^S,1:ndir)
55 double precision :: theta(ixI^S,1:ndir)
56 double precision :: rvertical(ixI^S),rhoadius(ixI^S)
57 double precision :: Itube
58 integer :: idir
59
60 itube=2.d0 &
61 *q_td99*l_td99*r_td99 &
62 /(l_td99**2+r_td99**2) &
63 /sqrt(l_td99**2+r_td99**2) &
64 /(log(8.d0*r_td99/a_td99)-1.5d0+li_td99*0.5d0)
65
66 rvertical(ixo^s)=dsqrt(x(ixo^s,2)**2+(x(ixo^s,3)+d_td99)**2)
67 rhoadius(ixo^s)=dsqrt(x(ixo^s,1)**2+(rvertical(ixo^s)-r_td99)**2)
68 theta(ixo^s,1)=0.d0
69 theta(ixo^s,2)=-(x(ixo^s,3)+d_td99)/rvertical(ixo^s)
70 theta(ixo^s,3)=x(ixo^s,2)/rvertical(ixo^s)
71
72 do idir=1,ndir
73 b(ixo^s,idir)= b(ixo^s,idir)+ &
74 izero_td99*2.0d0/sqrt(4.0d0*dpi)*theta(ixo^s,idir)*&
75 (1.0d0/sqrt(x(ixo^s,2)**2+(x(ixo^s,3)+d_td99)**2)-1.0d0/r_td99 &
76 +sqrt(1.d0/r_td99**2+2.0d0/a_td99**2*itube**2/izero_td99**2 &
77 *max(0.0d0,1.0d0-rhoadius(ixo^s)**2/a_td99**2)*p_bt_ratio))
78 end do
79
80 end subroutine td99_bth
81
82 subroutine td99_bi(ixI^L,ixO^L,x,B)
83 use mod_global_parameters
84 integer, intent(in) :: ixI^L, ixO^L
85 double precision, intent(in) :: x(ixI^S,1:ndim)
86 double precision, intent(inout) :: B(ixI^S,1:ndir)
87 double precision :: rvertical(ixI^S),rhoadius(ixI^S)
88
89 double precision :: Itube, Nt_TD99
90 double precision :: RperpUV(3),xUV(3)
91 double precision :: k,dkdx,dkdrv
92 double precision :: ka,dkadx,dkadrv
93 double precision :: Ak,dAkdk
94 double precision :: Aka,dAkdkka,d2Akdk2ka
95
96 double precision :: AIex,AIinwave,AI
97 double precision :: dAIexdx,dAIexdrv
98 double precision :: dAIinwavedx,dAIinwavedrv
99 double precision :: dAIdx,dAIdrv
100 integer :: i^D,idir
101
102 itube=2.d0 &
103 *q_td99*l_td99*r_td99 &
104 /(l_td99**2+r_td99**2) &
105 /sqrt(l_td99**2+r_td99**2) &
106 /(log(8.d0*r_td99/a_td99)-1.5d0+li_td99*0.5d0)
107
108 rvertical(ixo^s)=dsqrt(x(ixo^s,2)**2+(x(ixo^s,3)+d_td99)**2)
109 rhoadius(ixo^s)=dsqrt(x(ixo^s,1)**2+(rvertical(ixo^s)-r_td99)**2)
110
111 do i3 = ixomin3,ixomax3
112 do i2 = ixomin2,ixomax2
113 do i1 = ixomin1,ixomax1
114
115 rperpuv(1)=0.d0
116 rperpuv(2)=x(i^d,2)/rvertical(i^d)
117 rperpuv(3)=(x(i^d,3)+d_td99)/rvertical(i^d)
118 xuv(1)=1.d0
119 xuv(2)=0.d0
120 xuv(3)=0.d0
121
122 !Eq. 27 of TD99
123 k=2.d0*sqrt(rvertical(i^d)*r_td99 &
124 /((rvertical(i^d)+r_td99)**2+x(i^d,1)**2))
125 dkdx=-x(i^d,1)*k/((rvertical(i^d)+r_td99)**2+x(i^d,1)**2)
126 dkdrv=k*(r_td99**2+x(i^d,1)**2-rvertical(i^d)**2) &
127 /(2.d0*rvertical(i^d)*(r_td99**2+x(i^d,1)**2 &
128 +2.d0*r_td99*rvertical(i^d)+rvertical(i^d)**2))
129
130 !Eq. 30 of TD99
131 ka=2.d0*sqrt(rvertical(i^d)*r_td99 &
132 /(4.d0*rvertical(i^d)*r_td99+a_td99**2))
133 dkadx=0.d0
134 dkadrv=ka*a_td99**2 &
135 /(2.d0*rvertical(i^d)*(a_td99**2+4*r_td99*rvertical(i^d)))
136
137 !Eq. 26 of TD99
138 ak=1.d0/k*((2.d0-k*k)*ellf(dpi/2.d0,k) -2.d0*elle(dpi/2.d0,k))
139 !derivative of Eq. 26
140 dakdk=((2.d0-k*k)*elle(dpi/2.d0,k)-2.d0*(1.d0-k*k)*ellf(dpi/2.d0,k))/(k*k)/(1.d0-k*k)
141 !substitute ka into Eq. 26
142 aka=1.d0/ka*((2.d0-ka*ka)*ellf(dpi/2.d0,ka)-2.d0*elle(dpi/2.d0,ka))
143 !substitute ka into derivative of Eq. 26
144 dakdkka=((2.d0-ka*ka)*elle(dpi/2.d0,ka)- &
145 2.d0*(1.d0-ka*ka)*ellf(dpi/2.d0,ka))/(ka*ka)/(1.d0-ka*ka)
146 !substitute into derivative
147 d2akdk2ka=((7.d0*ka*ka-4.d0-ka**4)*elle(dpi/2.d0,ka)/(1.d0-ka*ka)+ &
148 (4.d0-5.d0*ka*ka)*ellf(dpi/2.d0,ka))/ka**3/(1.d0-ka*ka)
149
150 !Eq. 25 of TD99
151 aiex=itube*2.0d0*sqrt(r_td99/rvertical(i^d))*ak
152 !Eq. 28 of TD99
153 aiinwave=itube*2.0d0*sqrt(r_td99/rvertical(i^d))*(aka+dakdkka*(k-ka))
154 !Eq. 31 of TD99
155 if(a_td99 > rhoadius(i^d)) then
156 ai=aiinwave
157 else
158 ai=aiex
159 end if
160
161 !partial derivatives of Eq. 25 of TD99
162 daiexdx=itube*2.0d0*sqrt(r_td99/rvertical(i^d))*(dakdk*dkdx)
163 daiexdrv=itube*2.0d0*sqrt(r_td99/rvertical(i^d))*(dakdk*dkdrv)- &
164 aiex*0.5d0/rvertical(i^d)
165 !partial derivatives of Eq. 28 of TD99
166 daiinwavedx=itube*2.0d0*sqrt(r_td99/rvertical(i^d))*dakdkka*dkdx
167 daiinwavedrv=itube*2.0d0*sqrt(r_td99/rvertical(i^d)) &
168 *(dakdkka*dkdrv+d2akdk2ka*dkadrv*(k-ka))-aiinwave*0.5d0/rvertical(i^d)
169 !merge derivatives together
170 if(a_td99 > rhoadius(i^d)) then
171 daidx=daiinwavedx
172 daidrv=daiinwavedrv
173 else
174 daidx=daiexdx
175 daidrv=daiexdrv
176 end if
177
178 do idir=1,ndir
179 b(i^d,idir)=b(i^d,idir)+1.0d0/sqrt(4.0d0*dpi)* &
180 (-daidx*rperpuv(idir)+(daidrv+ai/rvertical(i^d))*xuv(idir))
181 end do
182
183 end do
184 end do
185 end do
186
187 end subroutine td99_bi
188}
189 ! numercial recipe's subroutines
190 function elle(phi,ak)
191 implicit none
192 double precision, intent(in) :: phi,ak
193 double precision :: elle
194 double precision :: cc,q,s
195 integer, parameter :: dp = kind(1.d0)
196 s=sin(phi)
197 cc=cos(phi)**2
198 q=(1.-s*ak)*(1.+s*ak)
199 elle=s*(rf(cc,q,1.0_dp)-((s*ak)**2)*rd(cc,q,1.0_dp)/3.0)
200 return
201 end function elle
202
203 function ellf(phi,ak)
204 double precision, intent(in) :: phi,ak
205 double precision :: ellf
206 double precision :: s
207 integer, parameter :: dp = kind(1.d0)
208 s=sin(phi)
209 ellf=s*rf(cos(phi)**2,(1.-s*ak)*(1.+s*ak),1.0_dp)
210 return
211 end function ellf
212
213 function ellpi(phi,en,ak)
214 double precision, intent(in) :: phi,en,ak
215 double precision :: ellpi
216 double precision :: cc,enss,q,s
217 integer, parameter :: dp = kind(1.d0)
218 s=sin(phi)
219 enss=en*s*s
220 cc=cos(phi)**2
221 q=(1.-s*ak)*(1.+s*ak)
222 ellpi=s*(rf(cc,q,1.0_dp)-enss*rj(cc,q,1.0_dp,1.0_dp+enss)/3.)
223 return
224 end function ellpi
225
226 FUNCTION rc(x,y)
227 double precision, intent(in) :: x,y
228 double precision :: rc
229 double precision :: errtol,tiny,sqrtny,big,tnbg,comp1,comp2,third, &
230 c1,c2,c3,c4
231 parameter(errtol=.04,tiny=1.69e-38,sqrtny=1.3e-19,big=3.e37, &
232 tnbg=tiny*big,comp1=2.236/sqrtny,comp2=tnbg*tnbg/25.,third=1./3., &
233 c1=.3,c2=1./7.,c3=.375,c4=9./22.)
234 double precision :: alamb,ave,s,w,xt,yt
235 if(x.lt.0..or.y.eq.0..or.(x+abs(y)).lt.tiny.or.(x+ &
236 abs(y)).gt.big.or.(y.lt.-comp1.and.x.gt.0..and.x.lt.comp2)) then
237 write(*,*)'invalid arguments in rc'
238 read(*,*)
239 endif
240 if(y.gt.0.)then
241 xt=x
242 yt=y
243 w=1.
244 else
245 xt=x-y
246 yt=-y
247 w=sqrt(x)/sqrt(xt)
248 endif
249 1 continue
250 alamb=2.*sqrt(xt)*sqrt(yt)+yt
251 xt=.25*(xt+alamb)
252 yt=.25*(yt+alamb)
253 ave=third*(xt+yt+yt)
254 s=(yt-ave)/ave
255 if(abs(s).gt.errtol)goto 1
256 rc=w*(1.+s*s*(c1+s*(c2+s*(c3+s*c4))))/sqrt(ave)
257 return
258 end function rc
259
260 FUNCTION rd(x,y,z)
261 double precision, intent(in) :: x,y,z
262 double precision :: rd
263 double precision :: errtol,tiny,big,c1,c2,c3,c4,c5,c6
264 parameter(errtol=.05,tiny=1.e-25,big=4.5e21,c1=3./14.,c2=1./6., &
265 c3=9./22.,c4=3./26.,c5=.25*c3,c6=1.5*c4)
266 double precision :: alamb,ave,delx,dely,delz,ea,eb,ec,ed,ee,fac,sqrtx,sqrty, &
267 sqrtz,sum,xt,yt,zt
268 if(min(x,y).lt.0..or.min(x+y,z).lt.tiny.or.max(x,y, &
269 z).gt.big) then
270 write(*,*)'invalid arguments in rd'
271 read(*,*)
272 endif
273 xt=x
274 yt=y
275 zt=z
276 sum=0.
277 fac=1.
278 1 continue
279 sqrtx=sqrt(xt)
280 sqrty=sqrt(yt)
281 sqrtz=sqrt(zt)
282 alamb=sqrtx*(sqrty+sqrtz)+sqrty*sqrtz
283 sum=sum+fac/(sqrtz*(zt+alamb))
284 fac=.25*fac
285 xt=.25*(xt+alamb)
286 yt=.25*(yt+alamb)
287 zt=.25*(zt+alamb)
288 ave=.2*(xt+yt+3.*zt)
289 delx=(ave-xt)/ave
290 dely=(ave-yt)/ave
291 delz=(ave-zt)/ave
292 if(max(abs(delx),abs(dely),abs(delz)).gt.errtol)goto 1
293 ea=delx*dely
294 eb=delz*delz
295 ec=ea-eb
296 ed=ea-6.*eb
297 ee=ed+ec+ec
298 rd=3.*sum+fac*(1.+ed*(-c1+c5*ed-c6*delz*ee)+delz*(c2*ee+delz*(-c3* &
299 ec+delz*c4*ea)))/(ave*sqrt(ave))
300 return
301 end function rd
302
303 FUNCTION rf(x,y,z)
304 double precision, intent(in) :: x,y,z
305 double precision :: rf
306 double precision :: errtol,tiny,big,third,c1,c2,c3,c4
307 parameter(errtol=.08,tiny=1.5e-38,big=3.e37,third=1./3., &
308 c1=1./24.,c2=.1,c3=3./44.,c4=1./14.)
309 double precision :: alamb,ave,delx,dely,delz,e2,e3,sqrtx,sqrty,sqrtz,xt,yt,zt
310 if(min(x,y,z).lt.0..or.min(x+y,x+z,y+z).lt.tiny.or.max(x,y, &
311 z).gt.big) then
312 write(*,*)'invalid arguments in rf'
313 read(*,*)
314 endif
315 xt=x
316 yt=y
317 zt=z
318 1 continue
319 sqrtx=sqrt(xt)
320 sqrty=sqrt(yt)
321 sqrtz=sqrt(zt)
322 alamb=sqrtx*(sqrty+sqrtz)+sqrty*sqrtz
323 xt=.25*(xt+alamb)
324 yt=.25*(yt+alamb)
325 zt=.25*(zt+alamb)
326 ave=third*(xt+yt+zt)
327 delx=(ave-xt)/ave
328 dely=(ave-yt)/ave
329 delz=(ave-zt)/ave
330 if(max(abs(delx),abs(dely),abs(delz)).gt.errtol)goto 1
331 e2=delx*dely-delz**2
332 e3=delx*dely*delz
333 rf=(1.+(c1*e2-c2-c3*e3)*e2+c4*e3)/sqrt(ave)
334 return
335 end function rf
336
337 FUNCTION rj(x,y,z,p)
338 double precision, intent(in) :: x,y,z,p
339 double precision :: rj
340 double precision :: errtol,tiny,big,c1,c2,c3,c4,c5,c6,c7,c8
341 parameter(errtol=.05,tiny=2.5e-13,big=9.e11,c1=3./14.,c2=1./3., &
342 c3=3./22.,c4=3./26.,c5=.75*c3,c6=1.5*c4,c7=.5*c2,c8=c3+c3)
343 double precision :: a,alamb,alpha,ave,b,beta,delp,delx,dely,delz,ea,eb,ec,ed,ee, &
344 fac,pt,rcx,rho,sqrtx,sqrty,sqrtz,sum,tau,xt,yt,zt
345 if(min(x,y,z).lt.0..or.min(x+y,x+z,y+z,abs(p)).lt.tiny.or.max(x,y, &
346 z,abs(p)).gt.big) then
347 write(*,*)'invalid arguments in rj'
348 read(*,*)
349 endif
350 sum=0.
351 fac=1.
352 if(p.gt.0.)then
353 xt=x
354 yt=y
355 zt=z
356 pt=p
357 else
358 xt=min(x,y,z)
359 zt=max(x,y,z)
360 yt=x+y+z-xt-zt
361 a=1./(yt-p)
362 b=a*(zt-yt)*(yt-xt)
363 pt=yt+b
364 rho=xt*zt/yt
365 tau=p*pt/yt
366 rcx=rc(rho,tau)
367 endif
368 1 continue
369 sqrtx=sqrt(xt)
370 sqrty=sqrt(yt)
371 sqrtz=sqrt(zt)
372 alamb=sqrtx*(sqrty+sqrtz)+sqrty*sqrtz
373 alpha=(pt*(sqrtx+sqrty+sqrtz)+sqrtx*sqrty*sqrtz)**2
374 beta=pt*(pt+alamb)**2
375 sum=sum+fac*rc(alpha,beta)
376 fac=.25*fac
377 xt=.25*(xt+alamb)
378 yt=.25*(yt+alamb)
379 zt=.25*(zt+alamb)
380 pt=.25*(pt+alamb)
381 ave=.2*(xt+yt+zt+pt+pt)
382 delx=(ave-xt)/ave
383 dely=(ave-yt)/ave
384 delz=(ave-zt)/ave
385 delp=(ave-pt)/ave
386 if(max(abs(delx),abs(dely),abs(delz),abs(delp)).gt.errtol)goto 1
387 ea=delx*(dely+delz)+dely*delz
388 eb=delx*dely*delz
389 ec=delp**2
390 ed=ea-3.*ec
391 ee=eb+2.*delp*(ea-ec)
392 rj=3.*sum+fac*(1.+ed*(-c1+c5*ed-c6*ee)+eb*(c7+delp*(-c8+delp*c4))+ &
393 delp*ea*(c2-delp*c3)-c2*delp*ec)/(ave*sqrt(ave))
394 if (p.le.0.) rj=a*(b*rj+3.*(rcx-rf(xt,yt,zt)))
395 return
396 end function rj
397end module mod_tdfluxrope
mod_global_parameters
This module contains definitions of global parameters and variables and some generic functions/subrou...
Definition mod_global_parameters.t:4
mod_global_parameters::d
double precision, dimension(:), allocatable, parameter d
Definition mod_global_parameters.t:23
mod_tdfluxrope
Definition mod_tdfluxrope.t:1
mod_tdfluxrope::l_td99
double precision l_td99
Definition mod_tdfluxrope.t:3
mod_tdfluxrope::d_td99
double precision d_td99
Definition mod_tdfluxrope.t:3
mod_tdfluxrope::q_td99
double precision q_td99
Definition mod_tdfluxrope.t:3
mod_tdfluxrope::izero_td99
double precision izero_td99
Definition mod_tdfluxrope.t:3
mod_tdfluxrope::rj
double precision function rj(x, y, z, p)
Definition mod_tdfluxrope.t:338
mod_tdfluxrope::elle
double precision function elle(phi, ak)
Definition mod_tdfluxrope.t:191
mod_tdfluxrope::td99
subroutine td99(ixil, ixol, x, b)
Definition mod_tdfluxrope.t:10
mod_tdfluxrope::ellf
double precision function ellf(phi, ak)
Definition mod_tdfluxrope.t:204
mod_tdfluxrope::rd
double precision function rd(x, y, z)
Definition mod_tdfluxrope.t:261
mod_tdfluxrope::td99_bi
subroutine td99_bi(ixil, ixol, x, b)
Definition mod_tdfluxrope.t:83
mod_tdfluxrope::td99_bth
subroutine td99_bth(ixil, ixol, x, b)
Definition mod_tdfluxrope.t:51
mod_tdfluxrope::td99_bq
subroutine td99_bq(ixil, ixol, x, b)
Definition mod_tdfluxrope.t:24
mod_tdfluxrope::li_td99
double precision li_td99
Definition mod_tdfluxrope.t:3
mod_tdfluxrope::rf
double precision function rf(x, y, z)
Definition mod_tdfluxrope.t:304
mod_tdfluxrope::r_td99
double precision r_td99
Definition mod_tdfluxrope.t:3
mod_tdfluxrope::p_bt_ratio
double precision p_bt_ratio
Definition mod_tdfluxrope.t:4
mod_tdfluxrope::a_td99
double precision a_td99
Definition mod_tdfluxrope.t:3
mod_tdfluxrope::ellpi
double precision function ellpi(phi, en, ak)
Definition mod_tdfluxrope.t:214
mod_tdfluxrope::rc
double precision function rc(x, y)
Definition mod_tdfluxrope.t:227
Generated by doxygen 1.9.8