Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

TLS for VectorPostprocessors #7819

Merged
merged 8 commits into from Oct 20, 2016
Merged

TLS for VectorPostprocessors #7819

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
59fd021
TLS for VectorPostprocessors
permcody Oct 5, 2016
2705db3
Revert "Fix thread fail (#7801)"
permcody Oct 7, 2016
2569c4d
Fix sign comparison warning
permcody Oct 7, 2016
119e35f
Rewrite SamplerBase to avoid temporaries
permcody Oct 7, 2016
2948396
Minor refactor of PointerSamplerBase
permcody Oct 12, 2016
1de886a
Fix broken test detected by assert
permcody Oct 12, 2016
7be7053
Update SphericalAverage VPP to avoid temporaries
permcody Oct 12, 2016
adcc4e9
Address comments
permcody Oct 13, 2016
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
2 changes: 1 addition & 1 deletion framework/include/base/FEProblem.h
View file
Open in desktop
Expand Up @@ -617,7 +617,7 @@ class FEProblem :
* Get the vectors for a specific VectorPostprocessor.
* @param vpp_name The name of the VectorPostprocessor
*/
const std::map<std::string, VectorPostprocessorValue*> & getVectorPostprocessorVectors(const std::string & vpp_name);
const std::map<std::string, VectorPostprocessorData::VectorPostprocessorState> & getVectorPostprocessorVectors(const std::string & vpp_name);

// Dampers /////
void addDamper(std::string damper_name, const std::string & name, InputParameters parameters);
Expand Down
14 changes: 5 additions & 9 deletions framework/include/vectorpostprocessors/PointSamplerBase.h
View file
Open in desktop
Expand Up @@ -45,10 +45,9 @@ class PointSamplerBase :
* Find the local element that contains the point. This will attempt to use a cached element to speed things up.
*
* @param p The point in physical space
* @param id A unique ID for this point.
* @return The Elem containing the point or NULL if this processor doesn't contain an element that contains this point.
*/
const Elem * getLocalElemContainingPoint(const Point & p, unsigned int id);
const Elem * getLocalElemContainingPoint(const Point & p);

/// The Mesh we're using
MooseMesh & _mesh;
Expand All @@ -59,17 +58,14 @@ class PointSamplerBase :
/// The ID to use for each point (yes, this is Real on purpose)
std::vector<Real> _ids;

/// Map of _points indices to the values
std::map<unsigned int, std::vector<Real> > _values;
/// Vector of values per point
std::vector<std::vector<Real> > _point_values;

/// Whether or not the Point was found on this processor (int because bool and uint don't work with MPI wrappers)
std::vector<int> _found_points;
/// Whether or not the Point was found on this processor (short because bool and char don't work with MPI wrappers)
std::vector<short> _found_points;

unsigned int _qp;

/// So we don't have to create and destroy this
std::vector<Point> _point_vec;

std::unique_ptr<PointLocatorBase> _pl;
};

Expand Down
14 changes: 1 addition & 13 deletions framework/include/vectorpostprocessors/SamplerBase.h
View file
Open in desktop
Expand Up @@ -94,7 +94,7 @@ class SamplerBase
std::vector<std::string> _variable_names;

/// What to sort by
unsigned int _sort_by;
const unsigned int _sort_by;

/// x coordinate of the points
VectorPostprocessorValue & _x;
Expand All @@ -107,18 +107,6 @@ class SamplerBase
VectorPostprocessorValue & _id;

std::vector<VectorPostprocessorValue *> _values;

/// x coordinate of the points
VectorPostprocessorValue _x_tmp;
/// y coordinate of the points
VectorPostprocessorValue _y_tmp;
/// x coordinate of the points
VectorPostprocessorValue _z_tmp;

/// The node ID of each point
VectorPostprocessorValue _id_tmp;

std::vector<VectorPostprocessorValue> _values_tmp;
};

#endif
5 changes: 1 addition & 4 deletions framework/include/vectorpostprocessors/SphericalAverage.h
View file
Open in desktop
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;
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Note to reviewer: I think we should avoid creating VectoPostprocessor instances in our objects. We may eventually want to make the constructor private like we have for Material and Variable containers to help prevent users from missing the reference in normal declarations. This particular vector is not a declared vector so I switched it to C++ vector.


/// aggregated global average vectors
std::vector<VectorPostprocessorValue *> _average;
Expand Down
5 changes: 5 additions & 0 deletions framework/include/vectorpostprocessors/VectorPostprocessor.h
View file
Open in desktop
Expand Up @@ -65,6 +65,11 @@ class VectorPostprocessor : public OutputInterface
FEProblem * _vpp_fe_problem;

friend class SamplerBase;

private:
THREAD_ID _vpp_tid;

std::map<std::string, VectorPostprocessorValue> _thread_local_vectors;
};

#endif
21 changes: 14 additions & 7 deletions framework/include/vectorpostprocessors/VectorPostprocessorData.h
View file
Open in desktop
Expand Up @@ -31,6 +31,12 @@ class VectorPostprocessorData : public Restartable
*/
VectorPostprocessorData(FEProblem & fe_problem);

struct VectorPostprocessorState
{
VectorPostprocessorValue * current;
VectorPostprocessorValue * old;
};

/**
* Initialization method, sets the current and old value to 0.0 for this
* VectorPostprocessor
Expand Down Expand Up @@ -64,26 +70,27 @@ class VectorPostprocessorData : public Restartable
/**
* Get the map of names -> VectorPostprocessor values. Exposed for error checking.
*/
const std::map<std::string, std::map<std::string, VectorPostprocessorValue*> > & values() const { return _values; }
// const std::map<std::string, std::map<std::string, VectorPostprocessorValue*> > & values() const { return _values; }

/**
* Get the map of vectors for a particular VectorPostprocessor
* @param vpp_name The name of the VectorPostprocessor
*/
const std::map<std::string, VectorPostprocessorValue*> & vectors(const std::string & vpp_name) { return _values[vpp_name]; }
const std::map<std::string, VectorPostprocessorState> & vectors(const std::string & vpp_name) const;

/**
* Copy the current post-processor values into old (i.e. shift it "back in time")
*/
void copyValuesBack();

protected:
private:
VectorPostprocessorValue & getVectorPostprocessorHelper(const VectorPostprocessorName & vpp_name, const std::string & vector_name, bool get_current);

/// Values of the post-processor at the current time
std::map<std::string, std::map<std::string, VectorPostprocessorValue*> > _values;
/// Values of the vector post-processor
std::map<std::string, std::map<std::string, VectorPostprocessorState> > _values;

/// Values of the post-processors at the time t-1
std::map<std::string, std::map<std::string, VectorPostprocessorValue*> > _values_old;
std::set<std::string> _requested_items;
std::set<std::string> _supplied_items;
};
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this for a planned dependency resolution?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes and no. Right now there's already dependency resolution on the objects themselves because they are UserObjects and sorted in the warehouse. The intention for this was to add in a new sanity check for the vectors themselves if you couple to vectors, there is no check on whether a vector that you couple to exists right now. I started on this but need to flesh it out quite a bit more and want to unify the way it works between VectorPostprocessors and Postprocessors. I can either continue to work on this or continue that work in another PR. I kind of want to handle it separately since this is already getting to be a substantial PR.


#endif //VECTORPOSTPROCESSORDATA_H
5 changes: 1 addition & 4 deletions framework/include/vectorpostprocessors/VolumeHistogram.h
View file
Open in desktop
Expand Up @@ -62,10 +62,7 @@ class VolumeHistogram : public ElementVectorPostprocessor
/// value mid point of the bin
VectorPostprocessorValue & _bin_center;

/// local thread copy of the volume vector
VectorPostprocessorValue _volume_tmp;

/// aggregated global volume vector
/// aggregated volume for the given bin
VectorPostprocessorValue & _volume;
};

Expand Down
2 changes: 1 addition & 1 deletion framework/src/base/FEProblem.C
View file
Open in desktop
Expand Up @@ -2311,7 +2311,7 @@ FEProblem::declareVectorPostprocessorVector(const VectorPostprocessorName & name
return _vpps_data.declareVector(name, vector_name);
}

const std::map<std::string, VectorPostprocessorValue*> &
const std::map<std::string, VectorPostprocessorData::VectorPostprocessorState> &
FEProblem::getVectorPostprocessorVectors(const std::string & vpp_name)
{
return _vpps_data.vectors(vpp_name);
Expand Down
12 changes: 8 additions & 4 deletions framework/src/outputs/TableOutput.C
View file
Open in desktop
Expand Up @@ -85,18 +85,22 @@ TableOutput::outputVectorPostprocessors()
// Loop through the postprocessor names and extract the values from the VectorPostprocessorData storage
for (const auto & vpp_name : out)
{
const std::map<std::string, VectorPostprocessorValue*> & vectors = _problem_ptr->getVectorPostprocessorVectors(vpp_name);
const auto & vectors = _problem_ptr->getVectorPostprocessorVectors(vpp_name);

FormattedTable & table = _vector_postprocessor_tables[vpp_name];
auto table_it = _vector_postprocessor_tables.lower_bound(vpp_name);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This set of 3 lines could use a bit of documentation. Took me a minute to figure out that you're just looking for or adding the table.

I mean: is this really worth removing the "accidental" insertion side effect? The old code was so simple...

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That's a fair criticism. I'm very, very tainted on RHS use of std::map::operator[]. I would say that I track down a bug about every 4-6 weeks that's related to careless use of that operator. I'd almost rather eradicate it's use anywhere in MOOSE but the syntax is very convenient which makes that hard to do.

This particular case is actually one of the few cases where RHS brackets are OK because we want to always create tables so I'm torn. In general, this is the syntax I'd like to move to because for the normal access case, you either add a straight up error when the item isn't found or often just an assertion when you are "sure" the item is there (yeah right...)

Here's another crazy idea... We use maps a lot in MOOSE. Maybe it's time to think about replacing our use of maps with a better container. We could start with just aliasing our container to std::map but quickly move to hiding implementation details such as the code you see there behind familiar operators like brackets. I've been using this lower_bound code lately because it's the best way to check for an item and insert if if necessary because you can use hint. Using find doesn't work because you always end up with a iterator to the end of the map meaning you have to do two tree searches to handle the insert case. Is it premature optimization? Maybe, but maybe not. I'm starting to put these in as I fix accidental insertion.

if (table_it == _vector_postprocessor_tables.end() || table_it->first != vpp_name)
table_it = _vector_postprocessor_tables.emplace_hint(table_it, vpp_name, FormattedTable());

FormattedTable & table = table_it->second;

table.clear();
table.outputTimeColumn(false);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yikes! Was this doing a copy before?!?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yep


for (const auto & vec_it : vectors)
{
VectorPostprocessorValue vector = *(vec_it.second);
const auto & vector = *vec_it.second.current;

for (unsigned int i=0; i<vector.size(); i++)
for (auto i = beginIndex(vector); i < vector.size(); ++i)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I like this!

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I do but I think we are in the minority. No other MOOSE devs are using this syntax, not even @dschwen. I wonder if we can do better still.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I used it in 513c5df after you made a comment to that effect :-D. It might just take a bit to catch on. I suggest we start converting loops (if we decide to stick with it). Seeing it in action in existing code may improve adoption.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This looks to me should be a method in vector, vector.begin_index(). This beginIndex looks ugly although correct.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah, sadly there is no method for creating an enumeration. There are other options and they are all uglier.

auto i = decltype(foo)(0)
or if you don't have a const type
decltype(foo) i = 0

table.addData(vec_it.first, vector[i], i);
}

Expand Down
42 changes: 25 additions & 17 deletions framework/src/vectorpostprocessors/PointSamplerBase.C
View file
Open in desktop
Expand Up @@ -35,8 +35,7 @@ PointSamplerBase::PointSamplerBase(const InputParameters & parameters) :
GeneralVectorPostprocessor(parameters),
CoupleableMooseVariableDependencyIntermediateInterface(this, false),
SamplerBase(parameters, this, _communicator),
_mesh(_subproblem.mesh()),
_point_vec(1) // Only going to evaluate one point at a time for now
_mesh(_subproblem.mesh())
{
std::vector<std::string> var_names(_coupled_moose_vars.size());

Expand All @@ -52,42 +51,51 @@ PointSamplerBase::initialize()
{
SamplerBase::initialize();

// We do this here just in case it's been destroyed and recreated becaue of mesh adaptivity.
// We do this here just in case it's been destroyed and recreated because of mesh adaptivity.
_pl = _mesh.getPointLocator();

// Reset the _found_points array
_found_points.resize(_points.size());
std::fill(_found_points.begin(), _found_points.end(), false);
// Reset the point arrays
_found_points.assign(_points.size(), false);

_point_values.resize(_points.size());
std::for_each(_point_values.begin(), _point_values.end(),
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Well done

[](std::vector<Real> & vec)
{
vec.clear();
});
}

void
PointSamplerBase::execute()
{
MeshTools::BoundingBox bbox = _mesh.getInflatedProcessorBoundingBox();

for (unsigned int i=0; i<_points.size(); i++)
/// So we don't have to create and destroy this
std::vector<Point> point_vec(1);

for (auto i = beginIndex(_points); i < _points.size(); ++i)
{
Point & p = _points[i];

// Do a bounding box check so we're not doing unnecessary PointLocator lookups
if (bbox.contains_point(p))
{
std::vector<Real> & values = _values[i];
auto & values = _point_values[i];

if (values.empty())
values.resize(_coupled_moose_vars.size());

// First find the element the hit lands in
const Elem * elem = getLocalElemContainingPoint(p, i);
const Elem * elem = getLocalElemContainingPoint(p);

if (elem)
{
// We have to pass a vector of points into reinitElemPhys
_point_vec[0] = p;
point_vec[0] = p;

_subproblem.reinitElemPhys(elem, _point_vec, 0); // Zero is for tid
_subproblem.reinitElemPhys(elem, point_vec, 0); // Zero is for tid

for (unsigned int j=0; j<_coupled_moose_vars.size(); j++)
for (auto j = beginIndex(_coupled_moose_vars); j < _coupled_moose_vars.size(); ++j)
values[j] = _coupled_moose_vars[j]->sln()[0]; // The zero is for the "qp"

_found_points[i] = true;
Expand All @@ -100,7 +108,7 @@ void
PointSamplerBase::finalize()
{
// Save off for speed
unsigned int pid = processor_id();
const auto pid = processor_id();
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Much better... we gotta get these unsigned ints out of the code... and, yes, I realize that a lot of them are my fault...


/*
* Figure out which processor is actually going "claim" each point.
Expand All @@ -111,27 +119,27 @@ PointSamplerBase::finalize()

_communicator.maxloc(_found_points, max_id);

for (unsigned int i=0; i<max_id.size(); i++)
for (auto i = beginIndex(max_id); i < max_id.size(); ++i)
{
// Only do this check on the proc zero because it's the same on every processor
// _found_points should contain all 1's at this point (ie every point was found by a proc)
if (pid == 0 && !_found_points[i])
mooseError("In " << name() << ", sample point not found: " << _points[i]);

if (max_id[i] == pid)
SamplerBase::addSample(_points[i], _ids[i], _values[i]);
SamplerBase::addSample(_points[i], _ids[i], _point_values[i]);
}

SamplerBase::finalize();
}

const Elem *
PointSamplerBase::getLocalElemContainingPoint(const Point & p, unsigned int /*id*/)
PointSamplerBase::getLocalElemContainingPoint(const Point & p)
{
const Elem * elem = (*_pl)(p);

if (elem && elem->processor_id() == processor_id())
return elem;

return NULL;
return nullptr;
}