1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
|
/*@@
@file GRHydro_AlfvenWaveM.F90
@date Oct 10, 2011
@author Bruno Mundim, Joshua Faber, Scott Noble
@desc
Circularly Polarized Alfven Wave test as implemented by
Beckwith and Stone Astrophys.J.Suppl. 193 (2011) 6, arXiv:1101.3573,
and Del Zanna et. al. A&A 473, 11 (2007), arXiv:0704.3206;
Other relevant references:
Stone et. al. Astrophys.J.Suppl. 178 (2008) 137, arXiv:0804.0402;
Gardiner and Stone JCP 227, 4123 (2008), arXiv:0712.2634;
Toth JCP 161, 605 (2000);
@enddesc
@@*/
#include "cctk.h"
#include "cctk_Parameters.h"
#include "cctk_Arguments.h"
#include "cctk_Functions.h"
#include "GRHydro_Macros.h"
#define velx(i,j,k) vel(i,j,k,1)
#define vely(i,j,k) vel(i,j,k,2)
#define velz(i,j,k) vel(i,j,k,3)
#define sx(i,j,k) scon(i,j,k,1)
#define sy(i,j,k) scon(i,j,k,2)
#define sz(i,j,k) scon(i,j,k,3)
#define Bvecx(i,j,k) Bvec(i,j,k,1)
#define Bvecy(i,j,k) Bvec(i,j,k,2)
#define Bvecz(i,j,k) Bvec(i,j,k,3)
#define Bconsx(i,j,k) Bcons(i,j,k,1)
#define Bconsy(i,j,k) Bcons(i,j,k,2)
#define Bconsz(i,j,k) Bcons(i,j,k,3)
/*@@
@routine GRHydro_AlfvenWaveM
@date Oct 10, 2011
@author Bruno Mundim, Joshua Faber, Scott Noble
@desc
Initial data for Circularly Polarized Alfven Wave test
@enddesc
@calls
@calledby
@history
Using GRHydro_AdvectedLoopM.F90 as a template.
@endhistory
@@*/
subroutine GRHydro_AlfvenWaveM(CCTK_ARGUMENTS)
implicit none
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_PARAMETERS
DECLARE_CCTK_FUNCTIONS
CCTK_INT :: i, j, k, nx, ny, nz
CCTK_REAL :: pi,gam,AA,wnbr
CCTK_REAL :: vxval,vyval,vzval, valf
CCTK_REAL :: rhoval,pressval,epsval,hval
CCTK_REAL :: Bxval, Byval, Bzval
CCTK_REAL :: dx,dy,dz
CCTK_REAL :: range_x,range_y,range_z,range_d
CCTK_REAL :: cos_theta, sin_theta
CCTK_REAL :: diaglength,xnew,vparallel,vperp,Bparallel,Bperp
CCTK_REAL :: Bvecx_d, Bvecz_d
CCTK_REAL :: velx_d, velz_d
CCTK_REAL :: det
CCTK_REAL :: t1,t2,t3
pi=4.0d0*atan(1.0d0)
!!$Adiabatic index for this test:
gam = (5.0d0/3.0d0)
!!$pressure, density, B^x, specific internal energy and enthalpy:
rhoval = 1.0d0
pressval = alfvenwave_pressure
Bxval = 1.0d0
epsval = pressval/(gam-1.0d0)/rhoval
hval = 1.0d0 + epsval + pressval/rhoval
!!$ DZ: rho=P=Bxval=AA = 1
!!$ Using DZ parameters, epsval=1.5, hval=3.5
!!$ Alfven Wave Amplitude:
AA=1.0d0
!!$ Alfven wave speed:
t1 = rhoval*hval+Bxval**2*(1.0d0+AA**2)
t2 = 2.0d0*AA*Bxval**2/t1
t3 = 0.5d0*(1.0d0+sqrt(1.0d0-t2**2))
valf = sqrt(Bxval**2/t1/t3)
write(*,*)'Alfven velocity:',valf
!!$ Using DZ parameters with P=1.0, we have:
!!$ t1=5.5
!!$ t2=4/11
!!$ t3=0.96577
!!$ valf=0.43389
!!$ Vx value:
vxval=0.0d0
nx = cctk_lsh(1)
ny = cctk_lsh(2)
nz = cctk_lsh(3)
dx = CCTK_DELTA_SPACE(1)
dy = CCTK_DELTA_SPACE(2)
dz = CCTK_DELTA_SPACE(3)
!!$ Note that the 3D test wasn't deviced to be used with AMR!
range_x = (cctk_gsh(1)-2*cctk_nghostzones(1))*dx
range_y = (cctk_gsh(2)-2*cctk_nghostzones(2))*dy
range_z = (cctk_gsh(3)-2*cctk_nghostzones(3))*dz
!!$ Alfven wave number
do i=1,nx
do j=1,ny
do k=1,nz
rho(i,j,k)=rhoval
press(i,j,k)=pressval
eps(i,j,k)=epsval
if (CCTK_EQUALS(alfvenwave_type,"1D")) then
wnbr = 2.0d0*pi/range_x
velx(i,j,k)=vxval
vely(i,j,k)=-valf*AA*cos(wnbr*x(i,j,k))
velz(i,j,k)=-valf*AA*sin(wnbr*x(i,j,k))
Bvecx(i,j,k)=Bxval
Bvecy(i,j,k)=AA*Bxval*cos(wnbr*x(i,j,k))
Bvecz(i,j,k)=AA*Bxval*sin(wnbr*x(i,j,k))
else if (CCTK_EQUALS(alfvenwave_type,"2D")) then
diaglength=range_x*range_y/range_d
range_d = sqrt(range_x**2+range_y**2)
cos_theta = range_y/range_d
sin_theta = range_x/range_d
wnbr = 2.0d0*pi/diaglength
xnew = cos_theta*x(i,j,k)+sin_theta*y(i,j,k)
vparallel=vxval
vperp=-valf*AA*cos(wnbr*xnew)
velx(i,j,k)=vparallel*cos_theta-vperp*sin_theta
vely(i,j,k)=vparallel*sin_theta+vperp*cos_theta
velz(i,j,k)=-valf*AA*sin(wnbr*xnew)
Bparallel=Bxval
Bperp=AA*Bxval*cos(wnbr*xnew)
Bvecx(i,j,k)=Bparallel*cos_theta-Bperp*sin_theta
Bvecy(i,j,k)=Bparallel*sin_theta+Bperp*cos_theta
Bvecz(i,j,k)=AA*Bxval*sin(wnbr*xnew)
else
call CCTK_WARN(0,"Alfven wave case not recognized!")
end if
det=SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k))
if (CCTK_EQUALS(GRHydro_eos_type,"Polytype")) then
call Prim2ConPolyM(GRHydro_eos_handle,gxx(i,j,k),gxy(i,j,k),&
gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),&
det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),&
tau(i,j,k),Bconsx(i,j,k),Bconsy(i,j,k),Bconsz(i,j,k),rho(i,j,k),&
velx(i,j,k),vely(i,j,k),velz(i,j,k),&
eps(i,j,k),press(i,j,k),Bvecx(i,j,k),Bvecy(i,j,k),Bvecz(i,j,k),&
w_lorentz(i,j,k))
else
call Prim2ConGenM(GRHydro_eos_handle,gxx(i,j,k),gxy(i,j,k),&
gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),&
det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),&
tau(i,j,k),Bconsx(i,j,k),Bconsy(i,j,k),Bconsz(i,j,k),rho(i,j,k),&
velx(i,j,k),vely(i,j,k),velz(i,j,k),&
eps(i,j,k),press(i,j,k),Bvecx(i,j,k),Bvecy(i,j,k),Bvecz(i,j,k),&
w_lorentz(i,j,k))
end if
enddo
enddo
enddo
densrhs = 0.d0
srhs = 0.d0
taurhs = 0.d0
Bconsrhs = 0.d0
return
end subroutine GRHydro_AlfvenWaveM
|
