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
|
#include <cctk.h>
#include <cctk_Arguments.h>
#include <cctk_Parameters.h>
#include "loopcontrol.h"
void Hydro_Analysis_PrepareReduction(CCTK_ARGUMENTS)
{
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_PARAMETERS
printf("local\n");
#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);
}
}
// Hydro_Analysis_rho_center_volume_weighted
void Hydro_Analysis_Reduction(CCTK_ARGUMENTS)
{
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_PARAMETERS
printf("global\n");
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_rho_max_search)
{
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 (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_WARN(0, "Error while reducing Hydro_Analysis_rho_center_volume_weighted_gf");
}
void Hydro_Analysis_LocationSearch(CCTK_ARGUMENTS)
{
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_PARAMETERS
printf("search\n");
/* Initialized MPI-local quantities */
CCTK_REAL local_rho_max_loc[3];
local_rho_max_loc[0] = 0.0;
local_rho_max_loc[1] = 0.0;
local_rho_max_loc[2] = 0.0;
/* 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];
}
}
}
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,
Hydro_Analysis_rho_max_loc, 3, CCTK_VARIABLE_REAL))
CCTK_WARN(0, "Error while reducing local_rho_max_loc");
}
|
