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
148
149
150
151
152
153
154
155
156
157
|
! Calculation of the sources for the level set function.
! $Header$
#include "cctk.h"
#include "cctk_Parameters.h"
#include "cctk_Arguments.h"
subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS)
use EHFinder_mod
implicit none
DECLARE_CCTK_PARAMETERS
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_FUNCTIONS
CCTK_INT :: i
CCTK_INT :: interp_handle, table_handle, status, coord_system_handle
CCTK_INT, dimension(1) :: lsh
CCTK_POINTER, dimension(3) :: interp_coords
CCTK_POINTER, dimension(14) :: out_arrays
CCTK_INT, dimension(12) :: in_arrays
CCTK_INT, dimension(14), parameter :: op_indices = (/ 0, 1, 2, 3, 4, &
5, 6, 7, 8, 9, &
10, 10, 10, 11 /), &
op_codes = (/ 0, 0, 0, 0, 0, &
0, 0, 0, 0, 0, &
1, 2, 3, 0 /)
CCTK_INT, dimension(14) :: out_types
CCTK_REAL :: alp2, dfux, dfuy, dfuz, factor
CCTK_REAL :: idetg, guxx, guxy, guxz, guyy, guyz, guzz
out_types = CCTK_VARIABLE_REAL
if ( CCTK_EQUALS ( metric_type, 'physical' ) ) then
call CCTK_GrouplshGN ( status, cctkGH, 1, lsh, "ehfinder::generators" )
if ( status .lt. 0 ) then
call CCTK_WARN ( 0, "cannot get local size for surface arrays" )
end if
! print*,lsh
call CCTK_InterpHandle ( interp_handle, &
"Lagrange polynomial interpolation" )
if ( interp_handle .lt. 0 ) then
call CCTK_WARN( 0, "Cannot get handle for interpolation. Forgot to activate an implementation providing interpolation operators??" )
end if
! For now order=2 is hard wired.
call Util_TableCreateFromString ( table_handle, "order=3" )
if ( table_handle .lt. 0 ) then
call CCTK_WARN( 0, "Cannot create parameter table for interpolator" )
end if
! Get the 3D coordinate system handle.
call CCTK_CoordSystemHandle ( coord_system_handle, "cart3d" )
if ( coord_system_handle .lt. 0) then
call CCTK_WARN( 0, "Cannot get handle for cart3d coordinate system. Forgot to activate an implementation providing coordinates ??" )
endif
interp_coords(1) = CCTK_PointerTo(xg)
interp_coords(2) = CCTK_PointerTo(yg)
interp_coords(3) = CCTK_PointerTo(zg)
out_arrays(1) = CCTK_PointerTo(alpg)
out_arrays(2) = CCTK_PointerTo(betaxg)
out_arrays(3) = CCTK_PointerTo(betayg)
out_arrays(4) = CCTK_PointerTo(betazg)
out_arrays(5) = CCTK_PointerTo(gxxg)
out_arrays(6) = CCTK_PointerTo(gxyg)
out_arrays(7) = CCTK_PointerTo(gxzg)
out_arrays(8) = CCTK_PointerTo(gyyg)
out_arrays(9) = CCTK_PointerTo(gyzg)
out_arrays(10) = CCTK_PointerTo(gzzg)
out_arrays(11) = CCTK_PointerTo(dfxg)
out_arrays(12) = CCTK_PointerTo(dfyg)
out_arrays(13) = CCTK_PointerTo(dfzg)
call CCTK_VarIndex ( in_arrays(1), "admbase::alp" )
call CCTK_VarIndex ( in_arrays(2), "admbase::betax" )
call CCTK_VarIndex ( in_arrays(3), "admbase::betay" )
call CCTK_VarIndex ( in_arrays(4), "admbase::betaz" )
call CCTK_VarIndex ( in_arrays(5), "admbase::gxx" )
call CCTK_VarIndex ( in_arrays(6), "admbase::gxy" )
call CCTK_VarIndex ( in_arrays(7), "admbase::gxz" )
call CCTK_VarIndex ( in_arrays(8), "admbase::gyy" )
call CCTK_VarIndex ( in_arrays(9), "admbase::gyz" )
call CCTK_VarIndex ( in_arrays(10), "admbase::gzz" )
call CCTK_VarIndex ( in_arrays(11), "ehfinder::f" )
call Util_TableSetIntArray ( status, table_handle, 13, &
op_indices(1:13), "operand_indices" )
if ( status .lt. 0 ) then
call CCTK_WARN ( 0, "Cannot set operand indices array in parameter table" )
endif
call Util_TableSetIntArray ( status, table_handle, 13, &
op_codes(1:13), "operation_codes" )
if ( status .lt. 0 ) then
call CCTK_WARN ( 0, "Cannot set operation codes array in parameter table" )
endif
call CCTK_InterpGridArrays ( status, cctkGH, 3, interp_handle, &
table_handle, coord_system_handle, &
lsh(1), CCTK_VARIABLE_REAL, &
interp_coords, 11, in_arrays(1:11), &
13, out_types(1:13), out_arrays(1:13) )
if ( status .lt. 0 ) then
call CCTK_INFO ( 'Interpolation failed.' )
end if
do i = 1, lsh(1)
alp2 = alpg(i)**2
guxx = gyyg(i) * gzzg(i) - gyzg(i)**2
guxy = gxzg(i) * gyzg(i) - gxyg(i) * gzzg(i)
guxz = gxyg(i) * gyzg(i) - gxzg(i) * gyyg(i)
idetg = one / ( gxxg(i) * guxx + gxyg(i) * guxy + gxzg(i) * guxz )
guxx = idetg * guxx
guxy = idetg * guxy
guxz = idetg * guxz
guyy = ( gxxg(i) * gzzg(i) - gxzg(i)**2 ) * idetg
guyz = ( gxyg(i) * gxzg(i) - gxxg(i) * gyzg(i) ) * idetg
guzz = ( gxxg(i) * gyyg(i) - gxyg(i)**2 ) * idetg
dfux = guxx * dfxg(i) + guxy * dfyg(i) + guxz * dfzg(i)
dfuy = guxy * dfxg(i) + guyy * dfyg(i) + guyz * dfzg(i)
dfuz = guxz * dfxg(i) + guyz * dfyg(i) + guzz * dfzg(i)
factor = alp2 / sqrt ( alp2 * ( dfux * dfxg(i) + &
dfuy * dfyg(i) + &
dfuz * dfzg(i) ) )
dxg(i) = - betaxg(i) + factor * dfux
dyg(i) = - betayg(i) + factor * dfuy
dzg(i) = - betazg(i) + factor * dfuz
! print*, alpg
! print*
! print*, betaxg, betayg, betazg
! print*
! print*, gxxg, gxyg, gxzg, gyyg, gyzg, gzzg
! print*
! print*, dfxg, dfyg, dfzg
! print*
! print*, alp2
! print*
! print*, dfux, dfuy, dfuz
! print*
! print*, factor
! print*
! print*, dxg(i), dyg(i), dzg(i)
! stop
end do
end if
return
end subroutine EHFinder_Generator_Sources
|
