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
|
/*@@
@file GRHydro_CylindricalExplosionM.F90
@date Apr 22, 2011
@author Scott Noble, Joshua Faber, Bruno Mundim
@desc
Cylindrical magnetized shocks.
@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_cylindricalexplosionM
@date Fri Apr 22 14:51:37 EDT 2011
@author Scott Noble, Joshua Faber, Bruno Mundim
@desc
Initial data for cylidrically symmetric shocks with
magnetic fields. Meant to recreate tests demonstrated
in Komissarov (1999).
@enddesc
@calls
@calledby
@history
Using GRHydro_ShockTubeM.F90 as a template.
@endhistory
@@*/
subroutine GRHydro_cylindricalexplosionM(CCTK_ARGUMENTS)
implicit none
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_PARAMETERS
DECLARE_CCTK_FUNCTIONS
CCTK_INT :: i, j, k, nx, ny, nz
CCTK_REAL :: direction, det
CCTK_REAL :: rhol, rhor, pressl, pressr
CCTK_REAL :: bvcxl,bvcyl,bvczl
CCTK_REAL :: tmp,tmp2,gam,r_inner,r_outer
CCTK_REAL :: cyl_fr
!!$ Check that user selected proper magnetic field ID. Warn if wrong one is
!!$ selected and proceed to ignore the value
if (.not. CCTK_EQUALS(initial_Bvec, "cylexp")) then
call CCTK_WARN(1, "When using the cyclindrical explosion initial data, please also select 'cylexp' for initial_Bvec. I will proceed as if this was set")
end if
!!$Uses Bx_init, By_init and Bz_init to set magnetic field strength
!!$Original tests had Bx_init = 0.1, 1.0 with By_init=Bz_init = 0
bvcxl = Bx_init
bvcyl = By_init
bvczl = Bz_init
!!$Inner radius and outer radius (Komissarov's defaults)
! r_inner = 8.d-1
! r_outer = 1.d0
r_inner = cyl_r_inner
r_outer = cyl_r_outer
!!$Inner values
rhor = cyl_rho_inner
pressr = cyl_press_inner
!!$Outer values
rhol = cyl_rho_outer
pressl = cyl_press_outer
!!$Adiabatic index for test
gam = gl_gamma
nx = cctk_lsh(1)
ny = cctk_lsh(2)
nz = cctk_lsh(3)
do i=1,nx
do j=1,ny
do k=1,nz
!!$direction represents the cylindrical radius here.
!!$TODO: maybe switch these over so that shocktube_tube choses the axis of the
!!$cylinder instead of the apparently random mapping
if (CCTK_EQUALS(shocktube_type,"xshock")) then
direction = sqrt((x(i,j,k)-shock_xpos)**2+&
(y(i,j,k)-shock_ypos)**2)
else if (CCTK_EQUALS(shocktube_type,"yshock")) then
direction = sqrt((y(i,j,k)-shock_ypos)**2+&
(z(i,j,k)-shock_zpos)**2)
else if (CCTK_EQUALS(shocktube_type,"zshock")) then
direction = sqrt((x(i,j,k)-shock_xpos)**2+&
(z(i,j,k)-shock_zpos)**2)
end if
tmp = cyl_fr(direction,r_inner,r_outer,rhol,rhor)
tmp2 = cyl_fr(direction,r_inner,r_outer,pressl,pressr)/( (gam - 1.d0) * tmp )
rho(i,j,k) = tmp
eps(i,j,k) = tmp2
velx(i,j,k) = 0.d0
vely(i,j,k) = 0.d0
velz(i,j,k) = 0.d0
Bvecx(i,j,k)=bvcxl
Bvecy(i,j,k)=bvcyl
Bvecz(i,j,k)=bvczl
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_CylindricalExplosionM
function cyl_fr(R,r_inner,r_outer,min_f,max_f)
implicit none
CCTK_REAL :: cyl_fr
CCTK_REAL :: R,r_inner,r_outer,min_f,max_f
if(R .gt. r_outer) then
cyl_fr = min_f
else if( R .lt. r_inner) then
cyl_fr = max_f
else
cyl_fr = exp( ( log(max_f)*(r_outer - R) + log(min_f)*(R - r_inner) ) / (r_outer - r_inner) )
end if
return
end function cyl_fr
|
