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
|
/*@@
@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 :: 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 = 1.0d0
Bxval = 1.0d0
epsval = pressval/(gam-1.0d0)/rhoval
hval = 1.0d0 + epsval + pressval/rhoval
!!$ 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)
!!$ 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)-1)*dx
range_y = (cctk_gsh(2)-1)*dy
range_z = (cctk_gsh(3)-1)*dz
range_d = sqrt(range_z**2+range_x**2)
cos_theta = range_z/range_d
sin_theta = range_x/range_d
!!$ Alfven wave number
wnbr = 2.0d0*pi/range_x
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
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))
!!$ if (CCTK_EQUALS(advectedloop_type,"3D")) then
!!$ Bvecx_d = cos_theta*Bvecx(i,j,k)+sin_theta*Bvecz(i,j,k)
!!$ Bvecz_d = cos_theta*Bvecz(i,j,k)-sin_theta*Bvecx(i,j,k)
!!$ velx_d = cos_theta*velx(i,j,k)+sin_theta*velz(i,j,k)
!!$ velz_d = cos_theta*velz(i,j,k)-sin_theta*velx(i,j,k)
!!$
!!$ Bvecx(i,j,k) = Bvecx_d
!!$ Bvecz(i,j,k) = Bvecz_d
!!$ velx(i,j,k) = velx_d
!!$ velz(i,j,k) = velz_d
!!$ 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
|
