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
|
/*@@
@file GRHydro_FluxM.F90
@date August 30, 2010
@author Joshua Faber, Scott Noble, Bruno Mundim, Pedro Montero, Ian Hawke
@desc
The routine to calculate the numerical flux function given a
specific state
@enddesc
@@*/
#include "cctk.h"
#include "cctk_Parameters.h"
#include "cctk_Arguments.h"
subroutine num_x_fluxM(dens,sx,sy,sz,tau,Bx,By,Bz,&
densf,sxf,syf,szf,tauf,Bxf,Byf,Bzf,vxt,vyt,vzt,pressstar,&
bsubx,bsuby,bsubz,ab0,w,det,alp,beta)
implicit none
CCTK_REAL :: dens,sx,sy,sz,tau,Bx,By,Bz
CCTK_REAL :: densf,sxf,syf,szf,tauf,Bxf,Byf,Bzf
CCTK_REAL :: vxt,vyt,vzt,bsubx,bsuby,bsubz,ab0,w
CCTK_REAL :: det,alp,beta,pressstar
CCTK_REAL :: velm
CCTK_REAL :: sdet,psipstar,psiBx,psiBy,psiBz
sdet=sqrt(det)
psipstar=pressstar*sdet
psiBx=Bx*sdet
psiBy=By*sdet
psiBz=Bz*sdet
!!$ We actually need all three values of vtilde = v^i - beta^i/alp, as well as
velm=vxt+beta/alp
!!$ GRHydro splits off alpha for later in the calculation, so we have psi^6 * [Anton eq.42]
!!$ In the notation of Anton et al.: [psi^6 D] vtilde^i
densf = dens * vxt
sxf = sx*vxt+psipstar-bsubx*psiBx/w
syf = sy*vxt-bsuby*psiBx/w
szf = sz*vxt-bsubz*psiBx/w
!!$ [psi^6 tau] vtilde^i +p* v^i - alp b^0 B^i/w
tauf = tau*vxt + psipstar*velm - ab0*psiBx/w
!!$ [psi^6 (B^k vtilde^i - B^i vtilde^k]
bxf = 0.0
byf = psiBy * vxt - psiBx*vyt
bzf = psiBz * vxt - psiBx*vzt
end subroutine num_x_fluxM
|
