aboutsummaryrefslogtreecommitdiff
path: root/src/GRHydro_Flux.F90
blob: 0b9eaab6d1443685f7c516adc91d873e2de31997 (plain)
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
/*@@
@file      GRHydro_Flux.F90
@date      Sat Jan 26 01:36:57 2002
@author    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"
 
 /*@@
   @routine    num_flux
   @date       Sat Jan 26 01:38:12 2002
   @author     Pedro Montero, Ian Hawke
   @desc 
   A general routine that works out the correct direction and returns
   the appropriate flux function
   @enddesc 
   @calls     
   @calledby   
   @history 
   Called from GR3D, original author Mark Miller.
   @endhistory 

@@*/
 
subroutine num_flux(xoffset,yoffset,zoffset,dens,sx,sy,sz,tau,&
     densf,sxf,syf,szf,tauf,velx,vely,velz,press,det,alp,beta)

  implicit none
  
  CCTK_REAL :: dens,sx,sy,sz,tau,velx,vely,velz,det
  CCTK_REAL :: densf,sxf,syf,szf,tauf
  CCTK_REAL :: alp,beta,press 
  integer :: xoffset,yoffset,zoffset

  densf = dens * ( (velx*xoffset + vely*yoffset + &
       velz*zoffset)  - beta / alp )
               
  sxf = sx * ( (velx*xoffset + vely*yoffset +&
       velz*zoffset)  - beta/alp) + &
       sqrt(det)*press*xoffset
               
  syf = sy * ( (velx*xoffset + vely*yoffset +&
       velz*zoffset)  - beta/alp) + &
       sqrt(det)*press*yoffset
              
  szf = sz * ( (velx*xoffset + vely*yoffset +&
       velz*zoffset)  - beta/alp ) + &
       sqrt(det)*press*zoffset
               
  tauf = tau * ( (velx*xoffset + vely*yoffset +&
       velz*zoffset)  - beta/alp ) + &
       sqrt(det)*press*(velx*xoffset +&
       vely*yoffset + velz*zoffset) 
 
end subroutine num_flux

 /*@@
   @routine    num_x_flux
   @date       Sat Jan 26 01:38:55 2002
   @author     Pedro Montero, Ian Hawke
   @desc 
   The numerical flux just along the x direction. Used for
   fluxes along the y,z direction by permuting the appropriate 
   arguments in the subroutine call.
   @enddesc 
   @calls     
   @calledby   
   @history 
   Culled and altered from GR3D, original author Mark Miller.
   @endhistory 

@@*/


subroutine num_x_flux(dens,sx,sy,sz,tau,&
     densf,sxf,syf,szf,tauf,vel,press,det,alp,beta)

  implicit none

  CCTK_REAL :: dens,sx,sy,sz,tau,vel,det
  CCTK_REAL :: densf,sxf,syf,szf,tauf
  CCTK_REAL :: alp,beta,press 
  CCTK_REAL :: velmbetainvalp,sqrtdetpress

  velmbetainvalp = vel - beta / alp
  sqrtdetpress = sqrt(det)*press
  densf = dens * ( velmbetainvalp )

  sxf = sx * velmbetainvalp + sqrtdetpress

  syf = sy * velmbetainvalp

  szf = sz * velmbetainvalp

  tauf = tau * velmbetainvalp + sqrtdetpress*vel
 
end subroutine num_x_flux