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
|
#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
*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 */
void Hydro_Analysis_LocationSearch(CCTK_ARGUMENTS)
{
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_PARAMETERS
/* 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
{
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. */
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 (fabs(level_rho_max_loc[3]) > 1.e-14) {
if (fabs(level_rho_max_loc[3]-1.) < 1.e-14)
{
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];
} else
CCTK_WARN(1, "Found more than one identical maximum, not setting anything.");
}
/* 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]); */
}
|
