MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
mod_tdfluxrope.t
Go to the documentation of this file.
2  implicit none
4  double precision :: p_bt_ratio
5 
6 contains
7 
8 {^ifthreed
9  subroutine td99(ixI^L,ixO^L,x,B)
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)
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)
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 &
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)
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
397 end module mod_tdfluxrope
This module contains definitions of global parameters and variables and some generic functions/subrou...
integer, dimension(:), allocatable, parameter d
subroutine td99_bq(ixIL, ixOL, x, B)
double precision l_td99
Definition: mod_tdfluxrope.t:3
double precision d_td99
Definition: mod_tdfluxrope.t:3
subroutine td99_bi(ixIL, ixOL, x, B)
double precision q_td99
Definition: mod_tdfluxrope.t:3
subroutine td99_bth(ixIL, ixOL, x, B)
double precision izero_td99
Definition: mod_tdfluxrope.t:3
double precision function rj(x, y, z, p)
double precision function elle(phi, ak)
subroutine td99(ixIL, ixOL, x, B)
double precision function ellf(phi, ak)
double precision function rd(x, y, z)
double precision li_td99
Definition: mod_tdfluxrope.t:3
double precision function rf(x, y, z)
double precision r_td99
Definition: mod_tdfluxrope.t:3
double precision p_bt_ratio
Definition: mod_tdfluxrope.t:4
double precision a_td99
Definition: mod_tdfluxrope.t:3
double precision function ellpi(phi, en, ak)
double precision function rc(x, y)