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
|
/*@@
@file GRHydro_PoloidalMagFieldM.F90
@date Oct 31, 2011
@author Bruno Mundim, Joshua Faber, Scott Noble
@desc
Poloidal Magnetic field implemented as in "General relativistic
simulations of magnetized binary neutron star mergers" -
by Yuk Tung Liu, Stuart L. Shapiro, Zachariah B. Etienne, and
Keisuke Taniguchi - Phys. Rev. D 78, 024012 (2008)
@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)
#define Avecx(i,j,k) Avec(i,j,k,1)
#define Avecy(i,j,k) Avec(i,j,k,2)
#define Avecz(i,j,k) Avec(i,j,k,3)
/*@@
@routine GRHydro_PoloidalMagFieldM
@date Oct 31, 2011
@author Bruno Mundim, Joshua Faber, Scott Noble
@desc
Poloidal Magnetic field implemented as in "General relativistic
simulations of magnetized binary neutron star mergers" -
by Yuk Tung Liu, Stuart L. Shapiro, Zachariah B. Etienne, and
Keisuke Taniguchi - Phys. Rev. D 78, 024012 (2008)
@enddesc
@calls
@calledby
@history
Using GRHydro_ShockTubeM.F90 as a template.
@endhistory
@@*/
subroutine GRHydro_PoloidalMagFieldM(CCTK_ARGUMENTS)
implicit none
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_PARAMETERS
DECLARE_CCTK_FUNCTIONS
CCTK_INT :: i, j, k, nx, ny, nz, set_Avec
CCTK_REAL :: det
CCTK_REAL :: sdet
CCTK_REAL :: dx,dy,dz
CCTK_REAL :: rhofac, delPcut, maxP_Pcut
CCTK_REAL :: AphiL, Ax, Ay, Az
CCTK_REAL :: rho_dx, rho_dy, rho_dz
CCTK_REAL :: press_dx, press_dy, press_dz
CCTK_REAL :: Aphi_dx, Aphi_dy, Aphi_dz
CCTK_REAL :: Ax_dy, Ax_dz, Ay_dx, Ay_dz
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)
set_Avec = 0
if ( CCTK_EQUALS(Bvec_evolution_method,"GRHydro_Avec") ) then
set_Avec = 1
end if
write(*,*)'GRHydro_InitData: Setting up initial poloidal magnetic field'
do i=2,nx-1
do j=2,ny-1
do k=2,nz-1
rhofac = 1.0d0-rho(i,j,k)/poloidal_rho_max
delPcut = press(i,j,k)-poloidal_P_cut
maxP_Pcut = max(delPcut,0.0d0)
AphiL = poloidal_A_b*rhofac**poloidal_n_p*maxP_Pcut
Ax = -y(i,j,k)*AphiL
Ay = x(i,j,k)*AphiL
Az = 0.0
!! write(*,*)'Before accessing rho(i,k,k)'
!! write(*,*)'rho(',i,',',j,',',k,') = ', rho(i,j,k)
!! write(*,*)'Ax, Ay, Az, Aphi = ', Ax, Ay, Az,AphiL
!! write(*,*)'rhofac = ', rhofac
!! write(*,*)'delPcut = ', delPcut
!! write(*,*)'maxP_Pcut = ', maxP_Pcut
rho_dx = 0.5d0*(rho(i+1,j,k)-rho(i-1,j,k))/dx
rho_dy = 0.5d0*(rho(i,j+1,k)-rho(i,j-1,k))/dy
rho_dz = 0.5d0*(rho(i,j,k+1)-rho(i,j,k-1))/dz
press_dx = 0.5d0*(press(i+1,j,k)-press(i-1,j,k))/dx
press_dy = 0.5d0*(press(i,j+1,k)-press(i,j-1,k))/dy
press_dz = 0.5d0*(press(i,j,k+1)-press(i,j,k-1))/dz
Aphi_dx = poloidal_A_b*rhofac**(poloidal_n_p-1)*maxP_Pcut* &
( rhofac*press_dx/delPcut - poloidal_n_p*rho_dx/poloidal_rho_max )
Aphi_dy = poloidal_A_b*rhofac**(poloidal_n_p-1)*maxP_Pcut* &
( rhofac*press_dy/delPcut - poloidal_n_p*rho_dy/poloidal_rho_max )
Aphi_dz = poloidal_A_b*rhofac**(poloidal_n_p-1)*maxP_Pcut* &
( rhofac*press_dz/delPcut - poloidal_n_p*rho_dz/poloidal_rho_max )
Ax_dy = -AphiL - y(i,j,k)*Aphi_dy
Ax_dz = -y(i,j,k)*Aphi_dz
Ay_dx = AphiL + x(i,j,k)*Aphi_dx
Ay_dz = x(i,j,k)*Aphi_dz
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))
sdet = sqrt(det)
Bvecx(i,j,k) = Ay_dz/sdet
Bvecy(i,j,k) = - Ax_dz/sdet
Bvecz(i,j,k) = (Ax_dy-Ay_dx)/sdet
if ( set_Avec.gt.0 ) then
Avecx(i,j,k) = Ax
Avecy(i,j,k) = Ay
Avecz(i,j,k) = Az
end if
!Bvecx(i,j,k) = 0.0d0
!Bvecy(i,j,k) = 0.0d0
!Bvecz(i,j,k) = 0.00000001/sdet
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
if (evolve_temper .ne. 0) then
call Prim2ConGenM_hot(GRHydro_eos_handle,GRHydro_reflevel,&
i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),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),temperature(i,j,k),y_e(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
end if
enddo
enddo
enddo
densrhs = 0.d0
srhs = 0.d0
taurhs = 0.d0
Bconsrhs = 0.d0
if (clean_divergence .ne. 0) then
psidcrhs =0.0
endif
return
end subroutine GRHydro_PoloidalMagFieldM
|
