Skip to content

Commit

Permalink
Update SphericalAverage VPP to avoid temporaries
Browse files Browse the repository at this point in the history
  • Loading branch information
permcody committed Oct 12, 2016
1 parent 22c6f93 commit 56cb41b
Show file tree
Hide file tree
Showing 2 changed files with 22 additions and 25 deletions.
5 changes: 1 addition & 4 deletions framework/include/vectorpostprocessors/SphericalAverage.h
Expand Up @@ -64,11 +64,8 @@ class SphericalAverage : public ElementVectorPostprocessor
/// value mid point of the bin
VectorPostprocessorValue & _bin_center;

/// local thread copy of the sum vectors
std::vector<VectorPostprocessorValue> _sum_tmp;

/// sample count per bin
VectorPostprocessorValue _count_tmp;
std::vector<unsigned int> _counts;

/// aggregated global average vectors
std::vector<VectorPostprocessorValue *> _average;
Expand Down
42 changes: 21 additions & 21 deletions framework/src/vectorpostprocessors/SphericalAverage.C
Expand Up @@ -37,35 +37,36 @@ SphericalAverage::SphericalAverage(const InputParameters & parameters) :
_values(_nvals),
_empty_bin_value(getParam<Real>("empty_bin_value")),
_bin_center(declareVector("radius")),
_sum_tmp(_nvals),
_count_tmp(_nbins),
_counts(_nbins),
_average(_nvals)
{
if (coupledComponents("variable") != 1)
mooseError("SphericalAverage works on exactly one coupled variable");

// Note: We associate the local variable "i" with nbins and "j" with nvals throughout.

// couple variables initialize vectors
for (unsigned int j = 0; j < _nvals; ++j)
for (auto j = beginIndex(_average); j < _nvals; ++j)
{
_values[j] = &coupledValue("variable", j);
_average[j] = &declareVector(getVar("variable", j)->name());
}

// initialize the bin center value vector
_bin_center.resize(_nbins);
for (unsigned i = 0; i < _nbins; ++i)
for (auto i = beginIndex(_counts); i < _nbins; ++i)
_bin_center[i] = (i + 0.5) * _deltaR;
}

void
SphericalAverage::initialize()
{
// reset the histogram
for (unsigned int j = 0; j < _nvals; ++j)
_sum_tmp[j].assign(_nbins, 0.0);
for (auto vec_ptr : _average)
vec_ptr->assign(_nbins, 0.0);

// reset bin counts
_count_tmp.assign(_nbins, 0.0);
_counts.assign(_nbins, 0.0);
}

void
Expand All @@ -75,31 +76,30 @@ SphericalAverage::execute()
for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
{
// compute target bin
int bin = computeDistance() / _deltaR;
auto bin = computeDistance() / _deltaR;

// add the volume contributed by the current quadrature point
if (bin >= 0 && bin < static_cast<int>(_nbins))
{
for (unsigned int j = 0; j < _nvals; ++j)
_sum_tmp[j][bin] += (*_values[j])[_qp];
for (auto j = decltype(_nvals)(0); j < _nvals; ++j)
(*_average[j])[bin] += (*_values[j])[_qp];

_count_tmp[bin]++;
_counts[bin]++;
}
}
}

void
SphericalAverage::finalize()
{
gatherSum(_count_tmp);
gatherSum(_counts);

for (unsigned int j = 0; j < _nvals; ++j)
for (auto j = beginIndex(_average); j < _nvals; ++j)
{
gatherSum(_sum_tmp[j]);
(*_average[j]).resize(_nbins);
gatherSum(*_average[j]);

for (unsigned int i = 0; i < _nbins; ++i)
(*_average[j])[i] = _count_tmp[i] > 0 ? _sum_tmp[j][i] / _count_tmp[i] : _empty_bin_value;
for (auto i = beginIndex(_counts); i < _nbins; ++i)
(*_average[j])[i] = _counts[i] > 0 ? (*_average[j])[i] / static_cast<Real>(_counts[i]) : _empty_bin_value;
}
}

Expand All @@ -108,12 +108,12 @@ SphericalAverage::threadJoin(const UserObject & y)
{
const SphericalAverage & uo = static_cast<const SphericalAverage &>(y);

for (unsigned int i = 0; i < _nbins; ++i)
for (auto i = beginIndex(_counts); i < _nbins; ++i)
{
_count_tmp[i] += uo._count_tmp[i];
_counts[i] += uo._counts[i];

for (unsigned int j = 0; j < _nvals; ++j)
_sum_tmp[j][i] += uo._sum_tmp[j][i];
for (auto j = beginIndex(_average); j < _nvals; ++j)
(*_average[j])[i] += (*uo._average[j])[i];
}
}

Expand Down

0 comments on commit 56cb41b

Please sign in to comment.