MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
mod_tdfluxrope.t
Go to the documentation of this file.
1 module 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 
6 contains
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
397 end 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::td99_bq
subroutine td99_bq(ixIL, ixOL, x, B)
Definition: mod_tdfluxrope.t:24
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::td99_bi
subroutine td99_bi(ixIL, ixOL, x, B)
Definition: mod_tdfluxrope.t:83
mod_tdfluxrope::q_td99
double precision q_td99
Definition: mod_tdfluxrope.t:3
mod_tdfluxrope::td99_bth
subroutine td99_bth(ixIL, ixOL, x, B)
Definition: mod_tdfluxrope.t:51
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::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.1