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
|
/*@@
@file GRHydro_CalcBcom.F90
@date Jan 01, 2012
@author Bruno Mundim
@histpry
Based on GRHydro_TmunuM.F90 file.
@desc
The calculation of the magnetic pressure and the comoving magnetic
field.
@enddesc
@@*/
#include "cctk.h"
#include "cctk_Arguments.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 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 bcomx(i,j,k) bcom(i,j,k,1)
#define bcomy(i,j,k) bcom(i,j,k,2)
#define bcomz(i,j,k) bcom(i,j,k,3)
subroutine GRHydro_CalcBcom(CCTK_ARGUMENTS)
implicit none
DECLARE_CCTK_ARGUMENTS
CCTK_REAL :: velxlow, velylow, velzlow
CCTK_REAL :: Bvecxlow,Bvecylow,Bveczlow
CCTK_REAL :: bdotv,b2,bxlow,bylow,bzlow,dum1,dum2
CCTK_INT :: i,j,k
CCTK_INT :: GRHydro_UseGeneralCoordinates
integer :: timelevels
call CCTK_ActiveTimeLevels(timelevels, cctkGH, "GRHydro::bcom")
if(timelevels.eq.0) then
call CCTK_WARN(0,"No storage for GRHydro::bcom")
end if
! Only compute if Tmunu is not computed! (Bcom is also computed with Tmunu)
if (stress_energy_state .ne. 0) then
return
endif
!$OMP PARALLEL DO PRIVATE(i,j,k,velxlow, velylow, velzlow,&
!$OMP Bvecxlow,Bvecylow,Bveczlow, bdotv,dum1,dum2,b2,bxlow,bylow,bzlow)
do k = 1, cctk_lsh(3)
do j = 1, cctk_lsh(2)
do i = 1, cctk_lsh(1)
call calc_vlow_blow(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),&
gyy(i,j,k),gyz(i,j,k),gzz(i,j,k), &
velx(i,j,k),vely(i,j,k),velz(i,j,k),Bvecx(i,j,k),Bvecy(i,j,k),Bvecz(i,j,k), &
velxlow,velylow,velzlow,Bvecxlow,Bvecylow,Bveczlow, &
bdotv,b2,dum1,dum2,bxlow,bylow,bzlow)
bcom_sq(i,j,k) = b2
bcom0(i,j,k) = w_lorentz(i,j,k)*bdotv/alp(i,j,k)
bcomx(i,j,k) = Bvecx(i,j,k)/w_lorentz(i,j,k) + bcom0(i,j,k)*(alp(i,j,k)*velx(i,j,k)-betax(i,j,k))
bcomy(i,j,k) = Bvecy(i,j,k)/w_lorentz(i,j,k) + bcom0(i,j,k)*(alp(i,j,k)*vely(i,j,k)-betay(i,j,k))
bcomz(i,j,k) = Bvecz(i,j,k)/w_lorentz(i,j,k) + bcom0(i,j,k)*(alp(i,j,k)*velz(i,j,k)-betaz(i,j,k))
end do
end do
end do
!$OMP END PARALLEL DO
return
end subroutine GRHydro_CalcBcom
|
