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
|
/*@@
@file GRHydro_BvecfromAvec
@date Aug 31, 2010
@author Tanja Bode
@desc
Calculate B^i (at cell center) from Avec (at edges)
@enddesc
@@*/
#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"
#include "cctk_Functions.h"
#include "GRHydro_Macros.h"
#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 Avecx(i,j,k) Avec(i,j,k,1)
#define Avecy(i,j,k) Avec(i,j,k,2)
#define Avecz(i,j,k) Avec(i,j,k,3)
! Second order f.d. for volume-centered quantity from volume-centered quantity
#define DIFF_X_2(q) (0.5d0 * (q(i+1,j,k) - q(i-1,j,k)) * idx)
#define DIFF_Y_2(q) (0.5d0 * (q(i,j+1,k) - q(i,j-1,k)) * idy)
#define DIFF_Z_2(q) (0.5d0 * (q(i,j,k+1) - q(i,j,k-1)) * idz)
! Fourth order f.d. for volume-centered quantity from volume-centered quantity
#define DIFF_X_4(q) ((-q(i+2,j,k) + 8.d0 * q(i+1,j,k) - 8.d0 * q(i-1,j,k) + \
q(i-2,j,k)) / 12.d0 * idx)
#define DIFF_Y_4(q) ((-q(i,j+2,k) + 8.d0 * q(i,j+1,k) - 8.d0 * q(i,j-1,k) + \
q(i,j-2,k)) / 12.d0 * idy)
#define DIFF_Z_4(q) ((-q(i,j,k+2) + 8.d0 * q(i,j,k+1) - 8.d0 * q(i,j,k-1) + \
q(i,j,k-2)) / 12.d0 * idz)
/*@@
@routine GRHydro_BvecfromAvec
@date Oct 24
@author Tanja Bode
@desc
Calculate B^i (at cell center) from Avec (at edges)
@enddesc
@calls
@calledby
@history
@endhistory
@@*/
subroutine GRHydro_BvecfromAvec(CCTK_ARGUMENTS)
implicit none
CCTK_REAL :: Az_y, Ay_z, Ax_z, Az_x, Ay_x, Ax_y
CCTK_REAL :: sdet, isdet
CCTK_REAL :: dx, dy, dz, idx, idy, idz
CCTK_INT :: i, j, k, nx, ny, nz, GRHydro_UseGeneralCoordinates
logical, allocatable, dimension (:,:,:) :: force_spatial_second_order
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_PARAMETERS
!!DECLARE_CCTK_FUNCTIONS
if (GRHydro_UseGeneralCoordinates(cctkGH).ne.0) then
call CCTK_WARN(0,"Bvec from Avec only in Cartesian at the moment.");
end if
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)
idx = 1.d0/dx
idy = 1.d0/dy
idz = 1.d0/dz
allocate (force_spatial_second_order(nx,ny,nz))
force_spatial_second_order = .FALSE.
if (spatial_order > 2) then
!$OMP PARALLEL DO PRIVATE(i, j, k)
do k = 1 + GRHydro_stencil, nz - GRHydro_stencil
do j = 1 + GRHydro_stencil, ny - GRHydro_stencil
do i = 1 + GRHydro_stencil, nx - GRHydro_stencil
if ((i < 3).or.(i > cctk_lsh(1) - 2).or. &
(j < 3).or.(j > cctk_lsh(2) - 2).or. &
(k < 3).or.(k > cctk_lsh(3) - 2) ) then
force_spatial_second_order(i,j,k) = .TRUE.
else if ( use_mask > 0 ) then
if (minval(emask(i-2:i+2,j-2:j+2,k-2:k+2)) < 0.75d0) then
force_spatial_second_order(i,j,k) = .TRUE.
end if
end if
end do
end do
end do
!$OMP END PARALLEL DO
end if
call CCTK_WARN(1,"Bvec from Avec start.");
!$OMP PARALLEL DO PRIVATE( sdet, isdet, local_spatial_order, &
!$OMP Az_y, Ay_z, Ax_z, Az_x, Ay_x, Ax_y)
do k = GRHydro_stencil, cctk_lsh(3)-GRHydro_stencil+1
do j = GRHydro_stencil, cctk_lsh(2)-GRHydro_stencil+1
do i = GRHydro_stencil, cctk_lsh(1)-GRHydro_stencil+1
local_spatial_order = spatial_order
if (force_spatial_second_order(i,j,k)) then
local_spatial_order = 2
end if
if (local_spatial_order .eq. 2) then
Az_y = DIFF_Y_2(Avecz)
Ay_z = DIFF_Z_2(Avecy)
Ax_z = DIFF_Z_2(Avecx)
Az_x = DIFF_X_2(Avecz)
Ay_x = DIFF_X_2(Avecy)
Ax_y = DIFF_Y_2(Avecx)
else
Az_y = DIFF_Y_4(Avecz)
Ay_z = DIFF_Z_4(Avecy)
Ax_z = DIFF_Z_4(Avecx)
Az_x = DIFF_X_4(Avecz)
Ay_x = DIFF_X_4(Avecy)
Ax_y = DIFF_Y_4(Avecx)
end if
sdet = sdetg(i,j,k)
isdet = 1.d0/sdet
Bvecx(i,j,k) = isdet*( Az_y - Ay_z )
Bvecy(i,j,k) = isdet*( Ax_z - Az_x )
Bvecz(i,j,k) = isdet*( Ay_x - Ax_y )
end do
end do
end do
!$OMP END PARALLEL DO
call CCTK_WARN(1,"Bvec from Avec end.");
end subroutine GRHydro_BvecfromAvec
|
