OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

changed: put ParallelIstlInformation in separate compile unit

Browse Source
This commit is contained in:
Arne Morten Kvarving 2022-08-17 10:37:00 +02:00
parent 6afb386070
commit 371b2592f5
9 changed files with 389 additions and 295 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
CMakeLists_files.cmake
View File

@ -52,6 +52,7 @@ list (APPEND MAIN_SOURCE_FILES
opm/simulators/linalg/FlexibleSolver5.cpp
opm/simulators/linalg/FlexibleSolver6.cpp
opm/simulators/linalg/MILU.cpp
opm/simulators/linalg/ParallelIstlInformation.cpp
opm/simulators/linalg/PropertyTree.cpp
opm/simulators/linalg/setupPropertyTree.cpp
opm/simulators/utils/PartiallySupportedFlowKeywords.cpp

1
opm/simulators/flow/BlackoilModelEbos.hpp
View File

@ -39,7 +39,6 @@
#include <opm/grid/UnstructuredGrid.h>
#include <opm/simulators/timestepping/SimulatorReport.hpp>
#include <opm/simulators/linalg/ParallelIstlInformation.hpp>
#include <opm/core/props/phaseUsageFromDeck.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/common/Exceptions.hpp>

3
opm/simulators/flow/countGlobalCells.cpp
View File

@ -21,9 +21,12 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#include "countGlobalCells.hpp"
#include <cassert>
#include <numeric>
namespace Opm {
namespace detail {

1
opm/simulators/flow/countGlobalCells.hpp
View File

@ -24,7 +24,6 @@
#ifndef OPM_COUNTGLOBALCELLS_HEADER_INCLUDED
#define OPM_COUNTGLOBALCELLS_HEADER_INCLUDED
#include <opm/simulators/linalg/ParallelIstlInformation.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <dune/grid/common/gridview.hh>

355
opm/simulators/linalg/ParallelIstlInformation.cpp Normal file
View File

@ -0,0 +1,355 @@
/*
Copyright 2014, 2015 Dr. Markus Blatt - HPC-Simulation-Software & Services
Copyright 2014, 2015 Statoil ASA
Copyright 2015 NTNU
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#if HAVE_MPI && HAVE_DUNE_ISTL
#include <opm/simulators/linalg/ParallelIstlInformation.hpp>
#include <dune/common/enumset.hh>
#include <opm/common/ErrorMacros.hpp>
#include <cstddef>
#include <exception>
#include <mpi.h>
#include <numeric>
namespace
{
template<class T>
class IndexSetInserter
{
public:
using ParallelIndexSet = T;
using LocalIndex = typename ParallelIndexSet::LocalIndex;
using GlobalIndex = typename ParallelIndexSet::GlobalIndex;
IndexSetInserter(ParallelIndexSet& indexSet, const GlobalIndex& component_size,
std::size_t local_component_size, std::size_t num_components)
: indexSet_(&indexSet), component_size_(component_size),
local_component_size_(local_component_size),
num_components_(num_components)
{}
void operator()(const typename ParallelIndexSet::IndexPair& pair)
{
for(std::size_t i = 0; i < num_components_; i++)
indexSet_->add(i * component_size_ + pair.global(),
LocalIndex(i * local_component_size_ + pair.local(),
pair.local().attribute()));
}
private:
ParallelIndexSet* indexSet_;
/// \brief The global number of unknowns per component/equation.
GlobalIndex component_size_;
/// \brief The local number of unknowns per component/equation.
std::size_t local_component_size_;
/// \brief The number of components/equations.
std::size_t num_components_;
};
/** \brief gather/scatter callback for communcation */
template<typename T>
struct CopyGatherScatter
{
using V = typename Dune::CommPolicy<T>::IndexedType;
static V gather(const T& a, std::size_t i)
{
return a[i];
}
static void scatter(T& a, V v, std::size_t i)
{
a[i] = v;
}
};
template<int I=0, typename... BinaryOperators, typename... ReturnValues>
typename std::enable_if<I == sizeof...(BinaryOperators), void>::type
computeGlobalReduction(const std::tuple<ReturnValues...>&,
std::tuple<BinaryOperators...>&,
std::tuple<ReturnValues...>&)
{}
template<int I=0, typename... BinaryOperators, typename... ReturnValues>
typename std::enable_if<I !=sizeof...(BinaryOperators), void>::type
computeGlobalReduction(const std::tuple<ReturnValues...>& receivedValues,
std::tuple<BinaryOperators...>& operators,
std::tuple<ReturnValues...>& values)
{
auto& val = std::get<I>(values);
val = std::get<I>(operators).localOperator()(val, std::get<I>(receivedValues));
computeGlobalReduction<I+1>(receivedValues, operators, values);
}
template<int I=0, typename... Containers, typename... BinaryOperators, typename... ReturnValues>
typename std::enable_if<I==sizeof...(Containers), void>::type
computeLocalReduction(const std::tuple<Containers...>&,
std::tuple<BinaryOperators...>&,
std::tuple<ReturnValues...>&,
const std::vector<double>&)
{}
template<int I=0, typename... Containers, typename... BinaryOperators, typename... ReturnValues>
typename std::enable_if<I!=sizeof...(Containers), void>::type
computeLocalReduction(const std::tuple<Containers...>& containers,
std::tuple<BinaryOperators...>& operators,
std::tuple<ReturnValues...>& values,
const std::vector<double>& ownerMask)
{
const auto& container = std::get<I>(containers);
if (container.size())
{
auto& reduceOperator = std::get<I>(operators);
// Eigen:Block does not support STL iterators!!!!
// Therefore we need to rely on the harder random-access
// property of the containers. But this should be save, too.
// Just commenting out code in the hope that Eigen might improve
// in this regard in the future.
//auto newVal = container.begin();
auto mask = ownerMask.begin();
auto& value = std::get<I>(values);
value = reduceOperator.getInitialValue();
for (auto endVal = ownerMask.end(); mask != endVal; /*++newVal,*/ ++mask )
{
value = reduceOperator(value, container[mask-ownerMask.begin()], *mask);
}
}
computeLocalReduction<I+1>(containers, operators, values, ownerMask);
}
}
namespace Opm
{
namespace
{
template<class T>
struct is_tuple
: std::integral_constant<bool, false>
{};
template<typename... T>
struct is_tuple<std::tuple<T...> >
: std::integral_constant<bool, true>
{};
}
ParallelISTLInformation::ParallelISTLInformation()
: indexSet_(new ParallelIndexSet),
remoteIndices_(new RemoteIndices(*indexSet_, *indexSet_, MPI_COMM_WORLD)),
communicator_(MPI_COMM_WORLD)
{}
ParallelISTLInformation::ParallelISTLInformation(MPI_Comm communicator)
: indexSet_(new ParallelIndexSet),
remoteIndices_(new RemoteIndices(*indexSet_, *indexSet_, communicator)),
communicator_(communicator)
{}
ParallelISTLInformation::
ParallelISTLInformation(const std::shared_ptr<ParallelIndexSet>& indexSet,
const std::shared_ptr<RemoteIndices>& remoteIndices,
MPI_Comm communicator)
: indexSet_(indexSet), remoteIndices_(remoteIndices), communicator_(communicator)
{}
ParallelISTLInformation::ParallelISTLInformation(const ParallelISTLInformation& other)
: indexSet_(other.indexSet_), remoteIndices_(other.remoteIndices_),
communicator_(other.communicator_)
{}
void ParallelISTLInformation::copyValuesTo(ParallelIndexSet& indexSet,
RemoteIndices& remoteIndices,
std::size_t local_component_size,
std::size_t num_components) const
{
ParallelIndexSet::GlobalIndex global_component_size = local_component_size;
if ( num_components > 1 )
{
ParallelIndexSet::GlobalIndex max_gi = 0;
// component the max global index
for( auto i = indexSet_->begin(), end = indexSet_->end(); i != end; ++i )
{
max_gi = std::max(max_gi, i->global());
}
global_component_size = max_gi+1;
global_component_size = communicator_.max(global_component_size);
}
indexSet.beginResize();
IndexSetInserter<ParallelIndexSet> inserter(indexSet, global_component_size,
local_component_size, num_components);
std::for_each(indexSet_->begin(), indexSet_->end(), inserter);
indexSet.endResize();
remoteIndices.rebuild<false>();
}
template<class T>
void ParallelISTLInformation::copyOwnerToAll(const T& source, T& dest) const
{
using AS = Dune::OwnerOverlapCopyAttributeSet;
using CopySet = Dune::EnumItem<AS, AS::copy>;
using OwnerSet = Dune::EnumItem<AS, AS::owner>;
using OverlapSet = Dune::EnumItem<AS, AS::overlap>;
using OwnerOverlapSet = Dune::Combine<OwnerSet, OverlapSet, AS::AttributeSet>;
using AllSet = Dune::Combine<OwnerOverlapSet, CopySet, AS::AttributeSet>;
OwnerSet sourceFlags;
AllSet destFlags;
Dune::Interface interface(communicator_);
if( !remoteIndices_->isSynced() )
{
remoteIndices_->rebuild<false>();
}
interface.build(*remoteIndices_,sourceFlags,destFlags);
Dune::BufferedCommunicator communicator;
communicator.template build<T>(interface);
communicator.template forward<CopyGatherScatter<T>>(source,dest);
communicator.free();
}
template<class T>
const std::vector<double>&
ParallelISTLInformation::updateOwnerMask(const T& container) const
{
if (!indexSet_)
{
OPM_THROW(std::runtime_error, "Trying to update owner mask without parallel information!");
}
if (static_cast<std::size_t>(container.size()) != ownerMask_.size())
{
ownerMask_.resize(container.size(), 1.);
for (const auto& i : *indexSet_)
{
if (i.local().attribute() != Dune::OwnerOverlapCopyAttributeSet::owner)
{
ownerMask_[i.local().local()] = 0.;
}
}
}
return ownerMask_;
}
template<typename Container, typename BinaryOperator, typename T>
void ParallelISTLInformation::computeReduction(const Container& container,
BinaryOperator binaryOperator,
T& value) const
{
if constexpr (is_tuple<Container>())
computeTupleReduction(container, binaryOperator, value);
else
{
std::tuple<const Container&> containers = std::tuple<const Container&>(container);
auto values = std::make_tuple(value);
auto operators = std::make_tuple(binaryOperator);
computeTupleReduction(containers, operators, values);
value = std::get<0>(values);
}
}
template<typename... Containers, typename... BinaryOperators, typename... ReturnValues>
void ParallelISTLInformation::computeTupleReduction(const std::tuple<Containers...>& containers,
std::tuple<BinaryOperators...>& operators,
std::tuple<ReturnValues...>& values) const
{
static_assert(std::tuple_size<std::tuple<Containers...> >::value ==
std::tuple_size<std::tuple<BinaryOperators...> >::value,
"We need the same number of containers and binary operators");
static_assert(std::tuple_size<std::tuple<Containers...> >::value ==
std::tuple_size<std::tuple<ReturnValues...> >::value,
"We need the same number of containers and return values");
if (std::tuple_size<std::tuple<Containers...> >::value == 0)
{
return;
}
// Copy the initial values.
std::tuple<ReturnValues...> init = values;
updateOwnerMask(std::get<0>(containers));
computeLocalReduction(containers, operators, values, ownerMask_);
std::vector<std::tuple<ReturnValues...> > receivedValues(communicator_.size());
communicator_.allgather(&values, 1, &(receivedValues[0]));
values = init;
for (auto& rval : receivedValues)
{
computeGlobalReduction(rval, operators, values);
}
}
template<class T1>
auto
accumulateMaskedValues(const T1& container, const std::vector<double>* maskContainer)
-> decltype(container[0]*(*maskContainer)[0])
{
decltype(container[0]*(*maskContainer)[0]) initial = 0;
if (maskContainer)
{
return std::inner_product(container.begin(), container.end(),
maskContainer->begin(), initial);
}
else
{
return std::accumulate(container.begin(), container.end(), initial);
}
}
template<class T> using C1 = std::vector<T>;
template<class T> using Ops1 = Reduction::MaskIDOperator<std::plus<T>>;
template<class T>
using C2 = std::tuple<std::vector<T>,
std::vector<T>,
std::vector<T>,
std::vector<T>,
std::vector<T>>;
template<class T>
using Ops2 = std::tuple<decltype(Reduction::makeGlobalSumFunctor<T>()),
decltype(Reduction::makeGlobalMaxFunctor<T>()),
decltype(Reduction::makeGlobalMinFunctor<T>()),
decltype(Reduction::makeInnerProductFunctor<T>()),
decltype(Reduction::makeLInfinityNormFunctor<T>())>;
template<class T>
using Vals2 = std::tuple<T,T,T,T,T>;
#define INSTANCE1(T) \
template void ParallelISTLInformation::computeReduction<C1<T>,Ops1<T>,T>(const C1<T>&,Ops1<T>,T&) const;
#define INSTANCE2(T) \
template void ParallelISTLInformation::computeReduction<C2<T>,Ops2<T>,Vals2<T>>(const C2<T>&,Ops2<T>,Vals2<T>&) const;
#define INSTANCE(T) \
INSTANCE1(T) \
INSTANCE2(T)
INSTANCE(int)
INSTANCE(float)
INSTANCE(std::size_t)
#endif
}

320
opm/simulators/linalg/ParallelIstlInformation.hpp
View File

@ -18,48 +18,25 @@
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_PARALLELISTLINFORMTION_HEADER_INCLUDED
#define OPM_PARALLELISTLINFORMTION_HEADER_INCLUDED
#ifndef OPM_PARALLELISTLINFORMATION_HEADER_INCLUDED
#define OPM_PARALLELISTLINFORMATION_HEADER_INCLUDED
#include <opm/grid/UnstructuredGrid.h>
#include <opm/common/ErrorMacros.hpp>
#include <any>
#include <exception>
#include <algorithm>
#include <functional>
#include <limits>
#include <numeric>
#include <type_traits>
#include <vector>
#if HAVE_MPI && HAVE_DUNE_ISTL
#include <opm/common/utility/platform_dependent/disable_warnings.h>
#include <mpi.h>
#include <algorithm>
#include <limits>
#include <memory>
#include <tuple>
#include <type_traits>
#include <dune/istl/owneroverlapcopy.hh>
#include <dune/common/parallel/interface.hh>
#include <dune/common/parallel/communicator.hh>
#include <dune/common/enumset.hh>
#include <opm/common/utility/platform_dependent/reenable_warnings.h>
#include <opm/simulators/utils/ParallelCommunication.hpp>
namespace Opm
{
namespace
{
template<class T>
struct is_tuple
: std::integral_constant<bool, false>
{};
template<typename... T>
struct is_tuple<std::tuple<T...> >
: std::integral_constant<bool, true>
{};
}
/// \brief Class that encapsulates the parallelization information needed by the
/// ISTL solvers.
@ -67,49 +44,42 @@ class ParallelISTLInformation
{
public:
/// \brief The type of the parallel index set used.
typedef Dune::OwnerOverlapCopyCommunication<int, int>::ParallelIndexSet ParallelIndexSet;
using ParallelIndexSet = Dune::OwnerOverlapCopyCommunication<int, int>::ParallelIndexSet;
/// \brief The type of the remote indices information used.
typedef Dune::OwnerOverlapCopyCommunication<int, int>::RemoteIndices RemoteIndices;
using RemoteIndices = Dune::OwnerOverlapCopyCommunication<int, int>::RemoteIndices;
/// \brief Constructs an empty parallel information object using MPI_COMM_WORLD
ParallelISTLInformation()
: indexSet_(new ParallelIndexSet),
remoteIndices_(new RemoteIndices(*indexSet_, *indexSet_, MPI_COMM_WORLD)),
communicator_(MPI_COMM_WORLD)
{}
ParallelISTLInformation();
/// \brief Constructs an empty parallel information object using a communicator.
/// \param communicator The communicator to use.
ParallelISTLInformation(MPI_Comm communicator)
: indexSet_(new ParallelIndexSet),
remoteIndices_(new RemoteIndices(*indexSet_, *indexSet_, communicator)),
communicator_(communicator)
{}
ParallelISTLInformation(MPI_Comm communicator);
/// \brief Constructs a parallel information object from the specified information.
/// \param indexSet The parallel index set to use.
/// \param remoteIndices The remote indices information to use.
/// \param communicator The communicator to use.
ParallelISTLInformation(const std::shared_ptr<ParallelIndexSet>& indexSet,
const std::shared_ptr<RemoteIndices>& remoteIndices,
MPI_Comm communicator)
: indexSet_(indexSet), remoteIndices_(remoteIndices), communicator_(communicator)
{}
MPI_Comm communicator);
/// \brief Copy constructor.
///
/// The information will be shared by the the two objects.
ParallelISTLInformation(const ParallelISTLInformation& other)
: indexSet_(other.indexSet_), remoteIndices_(other.remoteIndices_),
communicator_(other.communicator_)
{}
ParallelISTLInformation(const ParallelISTLInformation& other);
/// \brief Get a pointer to the underlying index set.
std::shared_ptr<ParallelIndexSet> indexSet() const
{
return indexSet_;
}
/// \brief Get a pointer to the remote indices information.
std::shared_ptr<RemoteIndices> remoteIndices() const
{
return remoteIndices_;
}
/// \brief Get the Collective MPI communicator that we use.
Parallel::Communication communicator() const
{
@ -119,69 +89,17 @@ public:
/// \param[out] indexSet The object to store the index set in.
/// \param[out] remoteIndices The object to store the remote indices information in.
void copyValuesTo(ParallelIndexSet& indexSet, RemoteIndices& remoteIndices,
std::size_t local_component_size = 0, std::size_t num_components = 1) const
{
ParallelIndexSet::GlobalIndex global_component_size = local_component_size;
if ( num_components > 1 )
{
ParallelIndexSet::GlobalIndex max_gi = 0;
// component the max global index
for( auto i = indexSet_->begin(), end = indexSet_->end(); i != end; ++i )
{
max_gi = std::max(max_gi, i->global());
}
global_component_size = max_gi+1;
global_component_size = communicator_.max(global_component_size);
}
indexSet.beginResize();
IndexSetInserter<ParallelIndexSet> inserter(indexSet, global_component_size,
local_component_size, num_components);
std::for_each(indexSet_->begin(), indexSet_->end(), inserter);
indexSet.endResize();
remoteIndices.rebuild<false>();
}
std::size_t local_component_size = 0,
std::size_t num_components = 1) const;
/// \brief Communcate the dofs owned by us to the other process.
///
/// Afterwards all associated dofs will contain the same data.
template<class T>
void copyOwnerToAll (const T& source, T& dest) const
{
typedef Dune::Combine<Dune::EnumItem<Dune::OwnerOverlapCopyAttributeSet::AttributeSet,Dune::OwnerOverlapCopyAttributeSet::owner>,Dune::EnumItem<Dune::OwnerOverlapCopyAttributeSet::AttributeSet,Dune::OwnerOverlapCopyAttributeSet::overlap>,Dune::OwnerOverlapCopyAttributeSet::AttributeSet> OwnerOverlapSet;
typedef Dune::EnumItem<Dune::OwnerOverlapCopyAttributeSet::AttributeSet,Dune::OwnerOverlapCopyAttributeSet::owner> OwnerSet;
typedef Dune::Combine<OwnerOverlapSet, Dune::EnumItem<Dune::OwnerOverlapCopyAttributeSet::AttributeSet,Dune::OwnerOverlapCopyAttributeSet::copy>,Dune::OwnerOverlapCopyAttributeSet::AttributeSet> AllSet;
OwnerSet sourceFlags;
AllSet destFlags;
Dune::Interface interface(communicator_);
if( !remoteIndices_->isSynced() )
{
remoteIndices_->rebuild<false>();
}
interface.build(*remoteIndices_,sourceFlags,destFlags);
Dune::BufferedCommunicator communicator;
communicator.template build<T>(interface);
communicator.template forward<CopyGatherScatter<T> >(source,dest);
communicator.free();
}
void copyOwnerToAll (const T& source, T& dest) const;
template<class T>
const std::vector<double>& updateOwnerMask(const T& container) const
{
if( ! indexSet_ )
{
OPM_THROW(std::runtime_error, "Trying to update owner mask without parallel information!");
}
if( static_cast<std::size_t>(container.size())!= ownerMask_.size() )
{
ownerMask_.resize(container.size(), 1.);
for( auto i=indexSet_->begin(), end=indexSet_->end(); i!=end; ++i )
{
if (i->local().attribute()!=Dune::OwnerOverlapCopyAttributeSet::owner)
{
ownerMask_[i->local().local()] = 0.;
}
}
}
return ownerMask_;
}
const std::vector<double>& updateOwnerMask(const T& container) const;
/// \brief Get the owner Mask.
///
@ -214,166 +132,14 @@ public:
/// \param value The initial value or a tuple of them.
template<typename Container, typename BinaryOperator, typename T>
void computeReduction(const Container& container, BinaryOperator binaryOperator,
T& value) const
{
computeReduction(container, binaryOperator, value, is_tuple<Container>());
}
T& value) const;
private:
/// \brief compute the reductions for tuples.
///
/// This is a helper function to prepare for calling computeTupleReduction.
template<typename Container, typename BinaryOperator, typename T>
void computeReduction(const Container& container, BinaryOperator binaryOperator,
T& value, std::integral_constant<bool,true>) const
{
computeTupleReduction(container, binaryOperator, value);
}
/// \brief compute the reductions for non-tuples.
///
/// This is a helper function to prepare for calling computeTupleReduction.
template<typename Container, typename BinaryOperator, typename T>
void computeReduction(const Container& container, BinaryOperator binaryOperator,
T& value, std::integral_constant<bool,false>) const
{
std::tuple<const Container&> containers=std::tuple<const Container&>(container);
auto values=std::make_tuple(value);
auto operators=std::make_tuple(binaryOperator);
computeTupleReduction(containers, operators, values);
value=std::get<0>(values);
}
/// \brief Compute the reductions for tuples.
template<typename... Containers, typename... BinaryOperators, typename... ReturnValues>
void computeTupleReduction(const std::tuple<Containers...>& containers,
std::tuple<BinaryOperators...>& operators,
std::tuple<ReturnValues...>& values) const
{
static_assert(std::tuple_size<std::tuple<Containers...> >::value==
std::tuple_size<std::tuple<BinaryOperators...> >::value,
"We need the same number of containers and binary operators");
static_assert(std::tuple_size<std::tuple<Containers...> >::value==
std::tuple_size<std::tuple<ReturnValues...> >::value,
"We need the same number of containers and return values");
if( std::tuple_size<std::tuple<Containers...> >::value==0 )
{
return;
}
// Copy the initial values.
std::tuple<ReturnValues...> init=values;
updateOwnerMask(std::get<0>(containers));
computeLocalReduction(containers, operators, values);
std::vector<std::tuple<ReturnValues...> > receivedValues(communicator_.size());
communicator_.allgather(&values, 1, &(receivedValues[0]));
values=init;
for( auto rvals=receivedValues.begin(), endvals=receivedValues.end(); rvals!=endvals;
++rvals )
{
computeGlobalReduction(*rvals, operators, values);
}
}
/// \brief TMP for computing the the global reduction after receiving the local ones.
///
/// End of recursion.
template<int I=0, typename... BinaryOperators, typename... ReturnValues>
typename std::enable_if<I == sizeof...(BinaryOperators), void>::type
computeGlobalReduction(const std::tuple<ReturnValues...>&,
std::tuple<BinaryOperators...>&,
std::tuple<ReturnValues...>&) const
{}
/// \brief TMP for computing the the global reduction after receiving the local ones.
template<int I=0, typename... BinaryOperators, typename... ReturnValues>
typename std::enable_if<I !=sizeof...(BinaryOperators), void>::type
computeGlobalReduction(const std::tuple<ReturnValues...>& receivedValues,
std::tuple<BinaryOperators...>& operators,
std::tuple<ReturnValues...>& values) const
{
auto& val=std::get<I>(values);
val = std::get<I>(operators).localOperator()(val, std::get<I>(receivedValues));
computeGlobalReduction<I+1>(receivedValues, operators, values);
}
/// \brief TMP for computing the the local reduction on the DOF that the process owns.
///
/// End of recursion.
template<int I=0, typename... Containers, typename... BinaryOperators, typename... ReturnValues>
typename std::enable_if<I==sizeof...(Containers), void>::type
computeLocalReduction(const std::tuple<Containers...>&,
std::tuple<BinaryOperators...>&,
std::tuple<ReturnValues...>&) const
{}
/// \brief TMP for computing the the local reduction on the DOF that the process owns.
template<int I=0, typename... Containers, typename... BinaryOperators, typename... ReturnValues>
typename std::enable_if<I!=sizeof...(Containers), void>::type
computeLocalReduction(const std::tuple<Containers...>& containers,
std::tuple<BinaryOperators...>& operators,
std::tuple<ReturnValues...>& values) const
{
const auto& container = std::get<I>(containers);
if( container.size() )
{
auto& reduceOperator = std::get<I>(operators);
// Eigen:Block does not support STL iterators!!!!
// Therefore we need to rely on the harder random-access
// property of the containers. But this should be save, too.
// Just commenting out code in the hope that Eigen might improve
// in this regard in the future.
//auto newVal = container.begin();
auto mask = ownerMask_.begin();
auto& value = std::get<I>(values);
value = reduceOperator.getInitialValue();
std::tuple<ReturnValues...>& values) const;
for( auto endVal=ownerMask_.end(); mask!=endVal;
/*++newVal,*/ ++mask )
{
value = reduceOperator(value, container[mask-ownerMask_.begin()], *mask);
}
}
computeLocalReduction<I+1>(containers, operators, values);
}
/** \brief gather/scatter callback for communcation */
template<typename T>
struct CopyGatherScatter
{
typedef typename Dune::CommPolicy<T>::IndexedType V;
static V gather(const T& a, std::size_t i)
{
return a[i];
}
static void scatter(T& a, V v, std::size_t i)
{
a[i] = v;
}
};
template<class T>
class IndexSetInserter
{
public:
typedef T ParallelIndexSet;
typedef typename ParallelIndexSet::LocalIndex LocalIndex;
typedef typename ParallelIndexSet::GlobalIndex GlobalIndex;
IndexSetInserter(ParallelIndexSet& indexSet, const GlobalIndex& component_size,
std::size_t local_component_size, std::size_t num_components)
: indexSet_(&indexSet), component_size_(component_size),
local_component_size_(local_component_size),
num_components_(num_components)
{}
void operator()(const typename ParallelIndexSet::IndexPair& pair)
{
for(std::size_t i = 0; i < num_components_; i++)
indexSet_->add(i * component_size_ + pair.global(),
LocalIndex(i * local_component_size_ + pair.local(),
pair.local().attribute()));
}
private:
ParallelIndexSet* indexSet_;
/// \brief The global number of unknowns per component/equation.
GlobalIndex component_size_;
/// \brief The local number of unknowns per component/equation.
std::size_t local_component_size_;
/// \brief The number of components/equations.
std::size_t num_components_;
};
std::shared_ptr<ParallelIndexSet> indexSet_;
std::shared_ptr<RemoteIndices> remoteIndices_;
Dune::CollectiveCommunication<MPI_Comm> communicator_;
@ -663,21 +429,6 @@ private:
namespace Opm
{
/// \brief Extracts the information about the data decomposition from the grid for dune-istl
///
/// In the case that grid is a parallel grid this method will query it to get the information
/// about the data decompoisition and convert it to the format expected by the linear algebra
/// of dune-istl.
/// \warn for UnstructuredGrid this function doesn't do anything.
/// \param anyComm The handle to store the information in. If grid is a parallel grid
/// then this will ecapsulate an instance of ParallelISTLInformation.
/// \param grid The grid to inspect.
inline void extractParallelGridInformationToISTL(std::any& anyComm, const UnstructuredGrid& grid)
{
(void)anyComm; (void)grid;
}
/// \brief Accumulates entries masked with 1.
/// \param container The container whose values to accumulate.
/// \param maskContainer null pointer or a pointer to a container
@ -687,19 +438,8 @@ inline void extractParallelGridInformationToISTL(std::any& anyComm, const Unstru
template<class T1>
auto
accumulateMaskedValues(const T1& container, const std::vector<double>* maskContainer)
-> decltype(container[0]*(*maskContainer)[0])
{
decltype(container[0]*(*maskContainer)[0]) initial = 0;
-> decltype(container[0]*(*maskContainer)[0]);
if( maskContainer )
{
return std::inner_product(container.begin(), container.end(), maskContainer->begin(),
initial);
}else
{
return std::accumulate(container.begin(), container.end(), initial);
}
}
} // end namespace Opm
#endif

1
opm/simulators/wells/RateConverter.hpp
View File

@ -24,7 +24,6 @@
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/grid/utility/RegionMapping.hpp>
#include <opm/simulators/linalg/ParallelIstlInformation.hpp>
#include <opm/simulators/wells/RegionAttributeHelpers.hpp>
#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
#include <dune/grid/common/gridenums.hh>

1
opm/simulators/wells/RegionAttributeHelpers.hpp
View File

@ -25,7 +25,6 @@
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/grid/utility/RegionMapping.hpp>
#include <opm/simulators/linalg/ParallelIstlInformation.hpp>
#include <dune/grid/common/gridenums.hh>
#include <algorithm>

1
opm/simulators/wells/RegionAverageCalculator.hpp
View File

@ -22,7 +22,6 @@
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/simulators/wells/RegionAttributeHelpers.hpp>
#include <opm/simulators/linalg/ParallelIstlInformation.hpp>
#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
#include <dune/grid/common/gridenums.hh>