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
|
/*@@
@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)
/*@@
@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
CCTK_REAL :: det
CCTK_REAL :: dx,dy,dz
CCTK_REAL :: rhofac, delPcut, maxP_Pcut
CCTK_REAL :: Aphi, 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)
do i=1,nx
do j=1,ny
do k=1,nz
rhofac = 1.0d0-rho(i,j,k)/poloidal_rho_max
delPcut = press(i,j,k)-poloidal_P_cut
maxP_Pcut = max(delPcut,0.0d0)
Aphi = poloidal_A_b*rhofac**poloidal_n_p*maxP_Pcut
Ax = -y(i,j,k)*Aphi
Ay = x(i,j,k)*Aphi
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,Aphi
!! 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 = -Aphi - y(i,j,k)*Aphi_dy
Ax_dz = y(i,j,k)*Aphi_dz
Ay_dx = Aphi + x(i,j,k)*Aphi_dx
Ay_dz = x(i,j,k)*Aphi_dz
Bvecx(i,j,k) = alp(i,j,k)*Ay_dz
Bvecy(i,j,k) = -alp(i,j,k)*Ax_dz
Bvecz(i,j,k) = alp(i,j,k)*(Ax_dy-Ay_dx)
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_PoloidalMagFieldM
|
