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
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
|
#include <math.h>
#include <cctk.h>
#include <cctk_Arguments.h>
#include <cctk_Parameters.h>
#include "loopcontrol.h"
void Hydro_Analysis_Init(CCTK_ARGUMENTS)
{
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_PARAMETERS
if (Hydro_Analysis_comp_rho_max || Hydro_Analysis_comp_rho_max_origin_distance)
*Hydro_Analysis_rho_max = -1.;
}
/* Populate grid functions to be reduced in Hydro_Analysis_Reduction */
void Hydro_Analysis_PrepareReduction(CCTK_ARGUMENTS)
{
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_PARAMETERS
if (Hydro_Analysis_comp_vol_weighted_center_of_mass)
{
#pragma omp parallel
{
LC_LOOP3(local_reduction, i,j,k,
cctk_nghostzones[0], cctk_nghostzones[1], cctk_nghostzones[2],
cctk_lsh[0] - cctk_nghostzones[0], cctk_lsh[1] - cctk_nghostzones[1],
cctk_lsh[2] - cctk_nghostzones[2], cctk_lsh[0], cctk_lsh[1], cctk_lsh[2])
{
CCTK_INT i3D = CCTK_GFINDEX3D(cctkGH, i, j, k);
CCTK_INT size= cctk_lsh[0] * cctk_lsh[1] * cctk_lsh[2];
Hydro_Analysis_rho_center_volume_weighted_gf[i3D+0*size] = rho[i3D] * x[i3D];
Hydro_Analysis_rho_center_volume_weighted_gf[i3D+1*size] = rho[i3D] * y[i3D];
Hydro_Analysis_rho_center_volume_weighted_gf[i3D+2*size] = rho[i3D] * z[i3D];
}
LC_ENDLOOP3(local_reduction);
}
}
}
/* Reduce grid functions */
void Hydro_Analysis_Reduction(CCTK_ARGUMENTS)
{
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_PARAMETERS
CCTK_INT handle_max = CCTK_ReductionHandle("maximum");
if (handle_max < 0)
CCTK_WARN(0, "Unable to get reduction handle 'max'.");
CCTK_INT handle_sum = CCTK_ReductionHandle("sum");
if (handle_sum < 0)
CCTK_WARN(0, "Unable to get reduction handle 'sum'.");
if (Hydro_Analysis_comp_rho_max)
{
if (CCTK_Reduce(cctkGH, -1, handle_max, 1, CCTK_VARIABLE_REAL,
Hydro_Analysis_rho_max, 1,
CCTK_VarIndex("HydroBase::rho")))
CCTK_WARN(0, "Error while reducing HydroBase::rho");
}
if (Hydro_Analysis_comp_vol_weighted_center_of_mass)
{
if (CCTK_Reduce(cctkGH, -1, handle_sum, 1, CCTK_VARIABLE_REAL,
Hydro_Analysis_rho_sum, 1,
CCTK_VarIndex("HydroBase::rho")))
CCTK_WARN(0, "Error while reducing HydroBase::rho");
if (CCTK_Reduce(cctkGH, -1, handle_sum, 1, CCTK_VARIABLE_REAL,
&(Hydro_Analysis_rho_center_volume_weighted[0]), 1,
CCTK_VarIndex(
"Hydro_Analysis::Hydro_Analysis_rho_center_volume_weighted_gf[0]")) ||
CCTK_Reduce(cctkGH, -1, handle_sum, 1, CCTK_VARIABLE_REAL,
&(Hydro_Analysis_rho_center_volume_weighted[1]), 1,
CCTK_VarIndex(
"Hydro_Analysis::Hydro_Analysis_rho_center_volume_weighted_gf[1]")) ||
CCTK_Reduce(cctkGH, -1, handle_sum, 1, CCTK_VARIABLE_REAL,
&(Hydro_Analysis_rho_center_volume_weighted[2]), 1,
CCTK_VarIndex(
"Hydro_Analysis::Hydro_Analysis_rho_center_volume_weighted_gf[2]")))
CCTK_WARN(0, "Error while reducing Hydro_Analysis_rho_center_volume_weighted_gf");
}
}
/* Look for the location of the density maximum */
/* RH: this is scheduled global loop-local but does MPI calls. It will hang if
* some processors have more components than others. */
void Hydro_Analysis_LocationSearch(CCTK_ARGUMENTS)
{
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_PARAMETERS
static int have_warned_about_multiple_maxima = 0;
int set_maximum_location, warned_about_multiple_local_maxima = 0;
/* Initialized MPI-local quantities */
CCTK_REAL local_rho_max_loc[4], level_rho_max_loc[4];
local_rho_max_loc[0] = 0.0;
local_rho_max_loc[1] = 0.0;
local_rho_max_loc[2] = 0.0;
local_rho_max_loc[3] = 0.0; // Information if it was found at all on this level
/* Look for the location of the global maximum.
* This algorithm will have problems when that occurs at more than
* one location. */
#pragma omp parallel
{
LC_LOOP3(location_reduction, i,j,k,
cctk_nghostzones[0], cctk_nghostzones[1], cctk_nghostzones[2],
cctk_lsh[0] - cctk_nghostzones[0], cctk_lsh[1] - cctk_nghostzones[1],
cctk_lsh[2] - cctk_nghostzones[2], cctk_lsh[0], cctk_lsh[1], cctk_lsh[2])
{
CCTK_INT i3D = CCTK_GFINDEX3D(cctkGH, i, j, k);
if (rho[i3D] == *Hydro_Analysis_rho_max &&
(!Hydro_Analysis_rho_max_loc_only_positive_x || x[i3D] >= 0.0))
{
#pragma omp critical
{
// have to already warn here when we still have the actual coordinates around
if (verbosity_level >= 1 && round(local_rho_max_loc[3]) >= 1.)
{
if(!have_warned_about_multiple_maxima && !warned_about_multiple_local_maxima)
CCTK_WARN(1, "Found more than one identical maximum in single patch.");
if (verbosity_level >= 2)
{
if (round(local_rho_max_loc[3]) == 1.) { // once we detect the second maximum, output the first as well
CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
"my candidate: (%g,%g,%g) with value %g.",
local_rho_max_loc[0],local_rho_max_loc[1],local_rho_max_loc[2],
*Hydro_Analysis_rho_max);
}
CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
"my candidate: (%g,%g,%g) with value %g.",
x[i3D], y[i3D], z[i3D],
*Hydro_Analysis_rho_max);
}
warned_about_multiple_local_maxima = 1;
}
local_rho_max_loc[0] += x[i3D];
local_rho_max_loc[1] += y[i3D];
local_rho_max_loc[2] += z[i3D];
local_rho_max_loc[3] += 1.;
}
}
}
LC_ENDLOOP3(location_reduction);
}
/* After this at least one MPI process should have the location set.
* The others should have a value of 0.0, so that the sum reduction
* gives the value on each MPI process. This however is not true sometimes
* and this is the point where the algorithm can fail. */
/* RH: this assumption is false if a processor owns more than a single patch
* */
CCTK_INT handle_sum = CCTK_LocalArrayReductionHandle("sum");
if (handle_sum < 0)
CCTK_WARN(0, "Unable to get reduction handle 'sum'.");
if (CCTK_ReduceLocArrayToArray1D(cctkGH, -1, handle_sum, &local_rho_max_loc,
level_rho_max_loc, 4, CCTK_VARIABLE_REAL))
CCTK_WARN(0, "Error while reducing local_rho_max_loc");
if (round(level_rho_max_loc[3]) > 0.) {
if (round(level_rho_max_loc[3]) == 1.)
{
set_maximum_location = 1;
} else {
if (Hydro_Analysis_average_multiple_maxima_locations)
{
level_rho_max_loc[0] /= level_rho_max_loc[3];
level_rho_max_loc[1] /= level_rho_max_loc[3];
level_rho_max_loc[2] /= level_rho_max_loc[3];
if (verbosity_level >= 1 && !have_warned_about_multiple_maxima)
{
CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
"Found more than one (%d) identical maximum, using average location (%g,%g,%g).",
(int)round(level_rho_max_loc[3]),
level_rho_max_loc[0], level_rho_max_loc[1],
level_rho_max_loc[2]);
have_warned_about_multiple_maxima = 1;
}
set_maximum_location = 1;
} else {
if (verbosity_level >= 1)
{
CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
"Found more than one (%d) identical maximum, not setting anything.",
(int)round(level_rho_max_loc[3]));
if (verbosity_level >= 2 && !warned_about_multiple_local_maxima)
{
CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
"my candidate: (%g,%g,%g) with value %g.",
local_rho_max_loc[0],local_rho_max_loc[1],local_rho_max_loc[2],
*Hydro_Analysis_rho_max);
}
}
set_maximum_location = 0;
}
}
} else {
set_maximum_location = 0;
}
if (set_maximum_location)
{
Hydro_Analysis_rho_max_loc[0] = level_rho_max_loc[0];
Hydro_Analysis_rho_max_loc[1] = level_rho_max_loc[1];
Hydro_Analysis_rho_max_loc[2] = level_rho_max_loc[2];
}
/* CCTK_VInfo(CCTK_THORNSTRING, "New location: %g,%g,%g",
Hydro_Analysis_rho_max_loc[0],
Hydro_Analysis_rho_max_loc[1],
Hydro_Analysis_rho_max_loc[2]); */
}
|
