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
158
159
|
#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Functions.h"
#include "cctk_Parameters.h"
subroutine qlm_import_surface (CCTK_ARGUMENTS, hn)
use cctk
use qlm_boundary
implicit none
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_FUNCTIONS
DECLARE_CCTK_PARAMETERS
integer :: hn
integer :: i, j, sn
character :: msg*1000
if (veryverbose/=0) then
call CCTK_INFO ("Importing surface shape")
end if
if ((surface_index(hn) < 0 .or. surface_index(hn) >= nsurfaces) .and. &
CCTK_EQUALS(surface_name(hn), "")) then
call CCTK_WARN (0, "Illegal spherical surface index specified")
end if
sn = sf_IdFromName(surface_index(hn), surface_name(hn)) + 1
if (verbose/=0 .or. veryverbose/=0) then
write (msg, '("Importing from spherical surface ",a," ",i4)') &
surface_name(hn), sn - 1
call CCTK_INFO (msg)
end if
if (qlm_calc_error(hn) == 0 .and. cctk_iteration > qlm_iteration(hn)) then
! Cycle time levels
qlm_have_valid_data_p_p(hn) = qlm_have_valid_data_p(hn)
qlm_have_valid_data_p (hn) = qlm_have_valid_data (hn)
qlm_have_killing_vector_p_p(hn) = qlm_have_killing_vector_p(hn)
qlm_have_killing_vector_p (hn) = qlm_have_killing_vector (hn)
qlm_iteration(hn) = cctk_iteration
qlm_time_p_p(hn) = qlm_time_p(hn)
qlm_time_p (hn) = qlm_time (hn)
qlm_radius_p_p(hn) = qlm_radius_p(hn)
qlm_radius_p (hn) = qlm_radius (hn)
qlm_origin_x_p_p(hn) = qlm_origin_x_p(hn)
qlm_origin_x_p (hn) = qlm_origin_x (hn)
qlm_origin_y_p_p(hn) = qlm_origin_y_p(hn)
qlm_origin_y_p (hn) = qlm_origin_y (hn)
qlm_origin_z_p_p(hn) = qlm_origin_z_p(hn)
qlm_origin_z_p (hn) = qlm_origin_z (hn)
qlm_shape_p_p(:,:,hn) = qlm_shape_p(:,:,hn)
qlm_shape_p (:,:,hn) = qlm_shape (:,:,hn)
end if
! Check for valid horizon data
if (sf_valid(sn) <= 0) then
if (verbose/=0) then
call CCTK_INFO ("No valid horizon data found")
end if
qlm_calc_error(hn) = 1
qlm_have_valid_data(hn) = 0
qlm_have_killing_vector(hn) = 0
return
end if
! Import surface description
qlm_calc_error(hn) = 0
qlm_have_valid_data(hn) = 0
qlm_have_killing_vector(hn) = 1
if (qlm_have_valid_data_p(hn) == 0) then
qlm_timederiv_order(hn) = 0
else if (qlm_have_valid_data_p_p(hn) == 0) then
qlm_timederiv_order(hn) = 1
else
qlm_timederiv_order(hn) = 2
end if
qlm_time(hn) = cctk_time
#warning "TODO: Ensure that the surface parameters don't change"
qlm_nghoststheta(hn) = sf_nghoststheta(sn)
qlm_nghostsphi (hn) = sf_nghostsphi(sn)
qlm_ntheta (hn) = sf_ntheta(sn)
qlm_nphi (hn) = sf_nphi(sn)
qlm_origin_x (hn) = sf_origin_x(sn)
qlm_origin_y (hn) = sf_origin_y(sn)
qlm_origin_z (hn) = sf_origin_z(sn)
qlm_origin_theta(hn) = sf_origin_theta(sn)
qlm_origin_phi (hn) = sf_origin_phi(sn)
qlm_delta_theta (hn) = sf_delta_theta(sn)
qlm_delta_phi (hn) = sf_delta_phi(sn)
if (veryverbose /= 0) then
write (msg, '("calc error : ",i6)') qlm_calc_error(hn)
call CCTK_INFO (msg)
write (msg, '("time deriv order: ",i6)') qlm_timederiv_order(hn)
call CCTK_INFO (msg)
write (msg, '("time : ",g16.6)') qlm_time(hn)
call CCTK_INFO (msg)
write (msg, '("nghosts : ",2i6)') qlm_nghoststheta(hn), qlm_nghostsphi(hn)
call CCTK_INFO (msg)
write (msg, '("n : ",2i6)') qlm_ntheta(hn), qlm_nphi(hn)
call CCTK_INFO (msg)
write (msg, '("origin : ",3g16.6)') &
qlm_origin_x(hn), qlm_origin_y(hn), qlm_origin_z(hn)
call CCTK_INFO (msg)
write (msg, '("origin : ",2g16.6)') &
qlm_origin_theta(hn), qlm_origin_phi(hn)
call CCTK_INFO (msg)
write (msg, '("delta : ",2g16.6)') &
qlm_delta_theta(hn), qlm_delta_phi(hn)
call CCTK_INFO (msg)
end if
if (qlm_ntheta(hn) > maxntheta .or. qlm_nphi(hn) > maxnphi) then
call CCTK_WARN (0, "Surface is too large")
end if
if (qlm_nghoststheta(hn)<1 .or. qlm_nghostsphi(hn)<1) then
call CCTK_WARN (0, "Not enough ghost zones for the horizon surface -- need at least 1")
end if
! Import the surface
! Calculate the coordinates
do j = 1, qlm_nphi(hn)
do i = 1, qlm_ntheta(hn)
qlm_shape(i,j,hn) = sf_radius(i,j,sn)
end do
end do
if (mod(int(qlm_ntheta(hn) - 2*qlm_nghoststheta(hn)),2) /= 1) then
! We need a grid point on the equator
call CCTK_WARN (0, "The number of interior grid points in theta direction must be odd")
end if
if (mod(int(qlm_nphi(hn) - 2*qlm_nghostsphi(hn)),4) /= 0) then
! We need grid points on the four major meridians
call CCTK_WARN (0, "The number of interior grid points in phi direction must be a multiple of four")
end if
end subroutine qlm_import_surface
|
