47 #ifndef OPM_CPGRIDDATA_HEADER
48 #define OPM_CPGRIDDATA_HEADER
51 #include <opm/grid/utility/platform_dependent/disable_warnings.h>
53 #include <dune/common/version.hh>
54 #include <dune/common/parallel/mpihelper.hh>
56 #include <dune/istl/owneroverlapcopy.hh>
59 #if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 7)
60 #include <dune/common/parallel/communication.hh>
62 #include <dune/common/parallel/collectivecommunication.hh>
64 #include <dune/common/parallel/indexset.hh>
65 #include <dune/common/parallel/interface.hh>
66 #include <dune/common/parallel/plocalindex.hh>
67 #if DUNE_VERSION_NEWER(DUNE_GRID, 2, 7)
68 #include <dune/common/parallel/variablesizecommunicator.hh>
70 #include <opm/grid/utility/VariableSizeCommunicator.hpp>
72 #include <dune/grid/common/gridenums.hh>
74 #include <opm/grid/utility/platform_dependent/reenable_warnings.h>
77 #include <opm/input/eclipse/EclipseState/Grid/NNC.hpp>
85 #include "OrientedEntityTable.hpp"
86 #include "DefaultGeometryPolicy.hpp"
89 #include <opm/grid/utility/OpmParserIncludes.hpp>
91 #include "Entity2IndexDataHandle.hpp"
92 #include "DataHandleWrappers.hpp"
93 #include "GlobalIdMapping.hpp"
107 class LevelGlobalIdSet;
108 class PartitionTypeIndicator;
109 template<
int,
int>
class Geometry;
110 template<
int>
class Entity;
111 template<
int>
class EntityRep;
115 template<
class T,
int i>
struct Mover;
133 #ifndef MAX_DATA_COMMUNICATED_PER_ENTITY
153 explicit CpGridData(MPIHelper::MPICommunicator comm);
160 int size(
int codim)
const;
163 int size (GeometryType type)
const
166 return size(3 - type.dim());
186 void readEclipseFormat(
const std::string& filename,
bool periodic_extension,
bool turn_normals =
false);
197 void processEclipseFormat(
const Opm::Deck& deck,
bool periodic_extension,
bool turn_normals =
false,
bool clip_z =
false,
198 const std::vector<double>& poreVolume = std::vector<double>());
211 std::vector<std::size_t>
processEclipseFormat(
const Opm::EclipseGrid* ecl_grid, Opm::EclipseState* ecl_state,
212 bool periodic_extension,
bool turn_normals =
false,
bool clip_z =
false,
bool pinchActive =
true);
225 std::array<std::set<std::pair<int, int>>, 2>& nnc,
226 bool remove_ij_boundary,
bool turn_normals,
bool pinchActive);
235 void getIJK(
int c, std::array<int,3>& ijk)
const
237 int gc = global_cell_[c];
238 ijk[0] = gc % logical_cartesian_size_[0]; gc /= logical_cartesian_size_[0];
239 ijk[1] = gc % logical_cartesian_size_[1];
240 ijk[2] = gc / logical_cartesian_size_[1];
244 void computeUniqueBoundaryIds();
251 return use_unique_boundary_ids_;
258 use_unique_boundary_ids_ = uids;
259 if (use_unique_boundary_ids_ && unique_boundary_ids_.empty()) {
260 computeUniqueBoundaryIds();
283 return logical_cartesian_size_;
291 const std::vector<int>& cell_part);
298 template<
class DataHandle>
299 void communicate(DataHandle& data, InterfaceType iftype, CommunicationDirection dir);
302 #if DUNE_VERSION_NEWER(DUNE_GRID, 2, 7)
304 using Communicator = VariableSizeCommunicator<>;
307 using Communicator = Opm::VariableSizeCommunicator<>;
311 using InterfaceMap = Communicator::InterfaceMap;
314 using CommunicationType = Dune::OwnerOverlapCopyCommunication<int,int>;
317 using ParallelIndexSet = CommunicationType::ParallelIndexSet;
320 using RemoteIndices = Dune::RemoteIndices<ParallelIndexSet>;
325 CommunicationType& cellCommunication()
333 const CommunicationType& cellCommunication()
const
338 ParallelIndexSet& cellIndexSet()
340 return cellCommunication().indexSet();
343 const ParallelIndexSet& cellIndexSet()
const
345 return cellCommunication().indexSet();
348 RemoteIndices& cellRemoteIndices()
350 return cellCommunication().remoteIndices();
353 const RemoteIndices& cellRemoteIndices()
const
355 return cellCommunication().remoteIndices();
359 #ifdef HAVE_DUNE_ISTL
361 typedef Dune::OwnerOverlapCopyAttributeSet::AttributeSet
AttributeSet;
370 return aquifer_cells_;
376 void populateGlobalCellIndexSet();
385 template<
class DataHandle>
386 void gatherData(DataHandle& data,
CpGridData* global_view,
396 template<
int codim,
class DataHandle>
397 void gatherCodimData(DataHandle& data,
CpGridData* global_data,
406 template<
class DataHandle>
407 void scatterData(DataHandle& data,
CpGridData* global_data,
408 CpGridData* distributed_data,
const InterfaceMap& cell_inf,
409 const InterfaceMap& point_inf);
418 template<
int codim,
class DataHandle>
419 void scatterCodimData(DataHandle& data,
CpGridData* global_data,
430 template<
int codim,
class DataHandle>
432 const Interface& interface);
442 template<
int codim,
class DataHandle>
444 const InterfaceMap& interface);
448 void computeGeometry(
CpGrid& grid,
450 const std::vector<int>& globalAquiferCells,
453 std::vector<int>& aquiferCells,
455 const std::vector< std::array<int,8> >& cell2Points);
471 std::vector< std::array<int,8> > cell_to_point_;
478 std::array<int, 3> logical_cartesian_size_;
485 std::vector<int> global_cell_;
491 typedef FieldVector<double, 3> PointType;
506 typedef MPIHelper::MPICommunicator MPICommunicator;
507 #if DUNE_VERSION_NEWER(DUNE_GRID, 2, 7)
508 using Communication = Dune::Communication<MPICommunicator>;
510 using CollectiveCommunication = Dune::CollectiveCommunication<MPICommunicator>;
511 using Communication = Dune::CollectiveCommunication<MPICommunicator>;
514 Communication ccobj_;
517 bool use_unique_boundary_ids_;
524 std::vector<double> zcorn;
527 std::vector<int> aquifer_cells_;
532 CommunicationType cell_comm_;
535 std::tuple<Interface,Interface,Interface,Interface,Interface> cell_interfaces_;
544 std::tuple<InterfaceMap,InterfaceMap,InterfaceMap,InterfaceMap,InterfaceMap>
557 template<
int>
friend class Entity;
558 template<
int>
friend class EntityRep;
559 friend class Intersection;
560 friend class PartitionTypeIndicator;
572 T& getInterface(InterfaceType iftype,
573 std::tuple<T,T,T,T,T>& interfaces)
578 return std::get<0>(interfaces);
580 return std::get<1>(interfaces);
582 return std::get<2>(interfaces);
584 return std::get<3>(interfaces);
586 return std::get<4>(interfaces);
588 OPM_THROW(std::runtime_error,
"Invalid Interface type was used during communication");
593 template<
int codim,
class DataHandle>
594 void CpGridData::communicateCodim(Entity2IndexDataHandle<DataHandle, codim>& data, CommunicationDirection dir,
595 const Interface& interface)
597 this->
template communicateCodim<codim>(data, dir, interface.interfaces());
600 template<
int codim,
class DataHandle>
601 void CpGridData::communicateCodim(Entity2IndexDataHandle<DataHandle, codim>& data_wrapper, CommunicationDirection dir,
602 const InterfaceMap& interface)
604 Communicator comm(ccobj_, interface);
606 if(dir==ForwardCommunication)
607 comm.forward(data_wrapper);
609 comm.backward(data_wrapper);
613 template<
class DataHandle>
615 CommunicationDirection dir)
618 if(data.contains(3,0))
621 communicateCodim<0>(data_wrapper, dir, getInterface(iftype, cell_interfaces_));
623 if(data.contains(3,3))
626 communicateCodim<3>(data_wrapper, dir, getInterface(iftype, point_interfaces_));
638 #include <opm/grid/cpgrid/Entity.hpp>
639 #include <opm/grid/cpgrid/Indexsets.hpp>
653 data=buffer_[index_++];
655 void write(
const T& data)
657 buffer_[index_++]=data;
663 void resize(std::size_t size)
665 buffer_.resize(size);
669 std::vector<T> buffer_;
670 typename std::vector<T>::size_type index_;
672 template<
class DataHandle,
int codim>
677 template<
class DataHandle>
680 BaseMover(DataHandle& data)
684 void moveData(
const E& from,
const E& to)
686 std::size_t size=data_.size(from);
688 data_.gather(buffer, from);
690 data_.scatter(buffer, to, size);
693 MoveBuffer<typename DataHandle::DataType> buffer;
697 template<
class DataHandle>
698 struct Mover<DataHandle,0> :
public BaseMover<DataHandle>
700 Mover<DataHandle,0>(DataHandle& data, CpGridData* gatherView,
701 CpGridData* scatterView)
702 : BaseMover<DataHandle>(data), gatherView_(gatherView), scatterView_(scatterView)
705 void operator()(std::size_t from_cell_index,std::size_t to_cell_index)
707 Entity<0> from_entity=Entity<0>(*gatherView_, from_cell_index,
true);
708 Entity<0> to_entity=Entity<0>(*scatterView_, to_cell_index,
true);
709 this->moveData(from_entity, to_entity);
711 CpGridData* gatherView_;
712 CpGridData* scatterView_;
715 template<
class DataHandle>
716 struct Mover<DataHandle,1> :
public BaseMover<DataHandle>
718 Mover<DataHandle,1>(DataHandle& data, CpGridData* gatherView,
719 CpGridData* scatterView)
720 : BaseMover<DataHandle>(data), gatherView_(gatherView), scatterView_(scatterView)
723 void operator()(std::size_t from_cell_index,std::size_t to_cell_index)
725 typedef typename OrientedEntityTable<0,1>::row_type row_type;
726 EntityRep<0> from_cell=EntityRep<0>(from_cell_index,
true);
727 EntityRep<0> to_cell=EntityRep<0>(to_cell_index,
true);
728 OrientedEntityTable<0,1>& table = gatherView_->cell_to_face_;
729 row_type from_faces=table.operator[](from_cell);
730 row_type to_faces=scatterView_->cell_to_face_[to_cell];
732 for(
int i=0; i<from_faces.size(); ++i)
733 this->moveData(from_faces[i], to_faces[i]);
735 CpGridData *gatherView_;
736 CpGridData *scatterView_;
739 template<
class DataHandle>
740 struct Mover<DataHandle,3> :
public BaseMover<DataHandle>
742 Mover<DataHandle,3>(DataHandle& data, CpGridData* gatherView,
743 CpGridData* scatterView)
744 : BaseMover<DataHandle>(data), gatherView_(gatherView), scatterView_(scatterView)
746 void operator()(std::size_t from_cell_index,std::size_t to_cell_index)
748 const std::array<int,8>& from_cell_points=
749 gatherView_->cell_to_point_[from_cell_index];
750 const std::array<int,8>& to_cell_points=
751 scatterView_->cell_to_point_[to_cell_index];
752 for(std::size_t i=0; i<8; ++i)
754 this->moveData(Entity<3>(*gatherView_, from_cell_points[i],
true),
755 Entity<3>(*scatterView_, to_cell_points[i],
true));
758 CpGridData* gatherView_;
759 CpGridData* scatterView_;
764 template<
class DataHandle>
765 void CpGridData::scatterData(DataHandle& data, CpGridData* global_data,
766 CpGridData* distributed_data,
const InterfaceMap& cell_inf,
767 const InterfaceMap& point_inf)
770 if(data.contains(3,0))
772 Entity2IndexDataHandle<DataHandle, 0> data_wrapper(*global_data, *distributed_data, data);
773 communicateCodim<0>(data_wrapper, ForwardCommunication, cell_inf);
775 if(data.contains(3,3))
777 Entity2IndexDataHandle<DataHandle, 3> data_wrapper(*global_data, *distributed_data, data);
778 communicateCodim<3>(data_wrapper, ForwardCommunication, point_inf);
783 template<
int codim,
class DataHandle>
784 void CpGridData::scatterCodimData(DataHandle& data, CpGridData* global_data,
785 CpGridData* distributed_data)
787 CpGridData *gather_view, *scatter_view;
788 gather_view=global_data;
789 scatter_view=distributed_data;
791 mover::Mover<DataHandle,codim> mover(data, gather_view, scatter_view);
794 for(
auto index=distributed_data->cellIndexSet().begin(),
795 end = distributed_data->cellIndexSet().end();
798 std::size_t from=index->global();
799 std::size_t to=index->local();
807 template<
int codim,
class T,
class F>
808 void visitInterior(CpGridData& distributed_data, T begin, T endit, F& func)
810 for(T it=begin; it!=endit; ++it)
812 Entity<codim> entity(distributed_data, it-begin,
true);
813 PartitionType pt = entity.partitionType();
814 if(pt==Dune::InteriorEntity)
821 template<
class DataHandle>
822 struct GlobalIndexSizeGatherer
824 GlobalIndexSizeGatherer(DataHandle& data_,
825 std::vector<int>& ownedGlobalIndices_,
826 std::vector<int>& ownedSizes_)
827 : data(data_), ownedGlobalIndices(ownedGlobalIndices_), ownedSizes(ownedSizes_)
830 template<
class T,
class E>
831 void operator()(T& i, E& entity)
833 ownedGlobalIndices.push_back(i);
834 ownedSizes.push_back(data.size(entity));
837 std::vector<int>& ownedGlobalIndices;
838 std::vector<int>& ownedSizes;
841 template<
class DataHandle>
844 DataGatherer(mover::MoveBuffer<typename DataHandle::DataType>& buffer_,
846 : buffer(buffer_), data(data_)
849 template<
class T,
class E>
850 void operator()(T& , E& entity)
852 data.gather(buffer, entity);
854 mover::MoveBuffer<typename DataHandle::DataType>& buffer;
860 template<
class DataHandle>
861 void CpGridData::gatherData(DataHandle& data, CpGridData* global_data,
862 CpGridData* distributed_data)
865 if(data.contains(3,0))
866 gatherCodimData<0>(data, global_data, distributed_data);
867 if(data.contains(3,3))
868 gatherCodimData<3>(data, global_data, distributed_data);
872 template<
int codim,
class DataHandle>
873 void CpGridData::gatherCodimData(DataHandle& data, CpGridData* global_data,
874 CpGridData* distributed_data)
878 const std::vector<int>& mapping =
879 distributed_data->global_id_set_->getMapping<codim>();
883 std::vector<int> owned_global_indices;
884 std::vector<int> owned_sizes;
885 owned_global_indices.reserve(mapping.size());
886 owned_sizes.reserve(mapping.size());
888 GlobalIndexSizeGatherer<DataHandle> gisg(data, owned_global_indices, owned_sizes);
889 visitInterior<codim>(*distributed_data, mapping.begin(), mapping.end(), gisg);
892 int no_indices=owned_sizes.size();
895 if ( owned_global_indices.empty() )
896 owned_global_indices.resize(1);
897 if ( owned_sizes.empty() )
898 owned_sizes.resize(1);
899 std::vector<int> no_indices_to_recv(distributed_data->ccobj_.size());
900 distributed_data->ccobj_.allgather(&no_indices, 1, &(no_indices_to_recv[0]));
904 std::vector<int> displ(distributed_data->ccobj_.size()+1, 0);
905 std::transform(displ.begin(), displ.end()-1, no_indices_to_recv.begin(), displ.begin()+1,
907 int global_size=displ[displ.size()-1];
908 std::vector<int> global_indices(global_size);
909 std::vector<int> global_sizes(global_size);
910 MPI_Allgatherv(&(owned_global_indices[0]), no_indices, MPITraits<int>::getType(),
911 &(global_indices[0]), &(no_indices_to_recv[0]), &(displ[0]),
912 MPITraits<int>::getType(),
913 distributed_data->ccobj_);
914 MPI_Allgatherv(&(owned_sizes[0]), no_indices, MPITraits<int>::getType(),
915 &(global_sizes[0]), &(no_indices_to_recv[0]), &(displ[0]),
916 MPITraits<int>::getType(),
917 distributed_data->ccobj_);
918 std::vector<int>().swap(owned_global_indices);
920 std::vector<int> no_data_send(distributed_data->ccobj_.size());
921 for(
typename std::vector<int>::iterator begin=no_data_send.begin(),
922 i=begin, end=no_data_send.end(); i!=end; ++i)
923 *i = std::accumulate(global_sizes.begin()+displ[i-begin],
924 global_sizes.begin()+displ[i-begin+1], std::size_t());
926 std::vector<int>().swap(owned_sizes);
929 std::transform(displ.begin(), displ.end()-1, no_data_send.begin(), displ.begin()+1,
930 std::plus<std::size_t>());
932 int no_data_recv = displ[displ.size()-1];
935 mover::MoveBuffer<typename DataHandle::DataType> local_data_buffer, global_data_buffer;
936 if ( no_data_send[distributed_data->ccobj_.rank()] )
938 local_data_buffer.resize(no_data_send[distributed_data->ccobj_.rank()]);
942 local_data_buffer.resize(1);
944 global_data_buffer.resize(no_data_recv);
946 DataGatherer<DataHandle> gatherer(local_data_buffer, data);
947 visitInterior<codim>(*distributed_data, mapping.begin(), mapping.end(), gatherer);
948 MPI_Allgatherv(&(local_data_buffer.buffer_[0]), no_data_send[distributed_data->ccobj_.rank()],
949 MPITraits<typename DataHandle::DataType>::getType(),
950 &(global_data_buffer.buffer_[0]), &(no_data_send[0]), &(displ[0]),
951 MPITraits<typename DataHandle::DataType>::getType(),
952 distributed_data->ccobj_);
953 Entity2IndexDataHandle<DataHandle, codim> edata(*global_data, data);
955 for(
int i=0; i< codim; ++i)
956 offset+=global_data->size(i);
958 typename std::vector<int>::const_iterator s=global_sizes.begin();
959 for(
typename std::vector<int>::const_iterator i=global_indices.begin(),
960 end=global_indices.end();
963 edata.scatter(global_data_buffer, *i-offset, *s);
[ provides Dune::Grid ]
Definition: CpGrid.hpp:210
Struct that hods all the data needed to represent a Cpgrid.
Definition: CpGridData.hpp:123
@ MAX_DATA_PER_CELL
The maximum data items allowed per cell (DUNE < 2.5.2)
Definition: CpGridData.hpp:140
int size(GeometryType type) const
number of leaf entities per geometry type in this process
Definition: CpGridData.hpp:163
void communicate(DataHandle &data, InterfaceType iftype, CommunicationDirection dir)
communicate objects for all codims on a given level
Definition: CpGridData.hpp:614
bool uniqueBoundaryIds() const
Is the grid currently using unique boundary ids?
Definition: CpGridData.hpp:249
CpGridData()
Constructor.
Definition: CpGridData.cpp:43
void readSintefLegacyFormat(const std::string &grid_prefix)
Read the Sintef legacy grid format ('topogeom').
Definition: readSintefLegacyFormat.cpp:66
int size(int codim) const
number of leaf entities per codim in this process
Definition: CpGridData.cpp:144
void readEclipseFormat(const std::string &filename, bool periodic_extension, bool turn_normals=false)
Read the Eclipse grid format ('grdecl').
~CpGridData()
Destructor.
Definition: CpGridData.cpp:93
void getIJK(int c, std::array< int, 3 > &ijk) const
Extract Cartesian index triplet (i,j,k) of an active cell.
Definition: CpGridData.hpp:235
const IndexSet & indexSet() const
Get the index set.
Definition: CpGridData.hpp:274
const std::array< int, 3 > & logicalCartesianSize() const
The logical cartesian size of the grid.
Definition: CpGridData.hpp:281
void distributeGlobalGrid(CpGrid &grid, const CpGridData &view_data, const std::vector< int > &cell_part)
Redistribute a global grid.
Definition: CpGridData.cpp:1553
void processEclipseFormat(const grdecl &input_data, Opm::EclipseState *ecl_state, std::array< std::set< std::pair< int, int >>, 2 > &nnc, bool remove_ij_boundary, bool turn_normals, bool pinchActive)
Read the Eclipse grid format ('grdecl').
Definition: processEclipseFormat.cpp:249
AttributeSet
The type of the set of the attributes.
Definition: CpGridData.hpp:364
void setUniqueBoundaryIds(bool uids)
Set whether we want to have unique boundary ids.
Definition: CpGridData.hpp:256
void writeSintefLegacyFormat(const std::string &grid_prefix) const
Write the Sintef legacy grid format ('topogeom').
Definition: writeSintefLegacyFormat.cpp:72
const std::vector< int > & sortedNumAquiferCells() const
Get sorted active cell indices of numerical aquifer.
Definition: CpGridData.hpp:368
const std::vector< double > & zcornData() const
Return the internalized zcorn copy from the grid processing, if no cells were adjusted during the min...
Definition: CpGridData.hpp:267
Definition: DefaultGeometryPolicy.hpp:49
const EntityVariable< cpgrid::Geometry< 3 - codim, 3 >, codim > & geomVector() const
Definition: DefaultGeometryPolicy.hpp:74
Wrapper that turns a data handle suitable for dune-grid into one based on integers instead of entitie...
Definition: Entity2IndexDataHandle.hpp:54
This class encapsulates geometry for both vertices, intersections and cells.
Definition: Geometry.hpp:71
The global id set for Dune.
Definition: Indexsets.hpp:325
Definition: Indexsets.hpp:198
Definition: Indexsets.hpp:54
Definition: Indexsets.hpp:257
Definition: PartitionTypeIndicator.hpp:50
Copyright 2019 Equinor AS.
Definition: CartesianIndexMapper.hpp:10
Holds the implementation of the CpGrid as a pimple.
Definition: CellQuadrature.hpp:29
Low-level corner-point processing routines and supporting data structures.
Definition: CpGridData.hpp:115
Raw corner-point specification of a particular geological model.
Definition: preprocess.h:56