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
|
/*@@
@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
sdet=sqrt(det)
psipstar=pressstar*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*Bx/w
syf = sy*vxt-bsuby*Bx/w
szf = sz*vxt-bsubz*Bx/w
!!$ [psi^6 tau] vtilde^i +p* v^i - alp b^0 B^i/w
tauf = tau*vxt + psipstar*velm - ab0*Bx/w
!!$ [psi^6 (B^k vtilde^i - B^i vtilde^k)]
bxf = 0.0
byf = By * vxt - Bx*vyt
bzf = Bz * vxt - Bx*vzt
end subroutine num_x_fluxM
|
