MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_CoordinatesTransferFactory_kokkos_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // MueLu: A package for multigrid based preconditioning
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact
39 // Jonathan Hu (jhu@sandia.gov)
40 // Andrey Prokopenko (aprokop@sandia.gov)
41 // Ray Tuminaro (rstumin@sandia.gov)
42 //
43 // ***********************************************************************
44 //
45 // @HEADER
46 #ifndef MUELU_COORDINATESTRANSFER_FACTORY_KOKKOS_DEF_HPP
47 #define MUELU_COORDINATESTRANSFER_FACTORY_KOKKOS_DEF_HPP
48 
50 
51 #include <Xpetra_ImportFactory.hpp>
52 #include <Xpetra_IO.hpp>
53 #include <Xpetra_MultiVectorFactory.hpp>
54 #include <Xpetra_MapFactory.hpp>
55 
56 #include "MueLu_Aggregates_kokkos.hpp"
57 #include "MueLu_Level.hpp"
58 #include "MueLu_Monitor.hpp"
59 #include "MueLu_Utilities_kokkos.hpp"
60 
61 namespace MueLu {
62 
63  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class DeviceType>
65  RCP<ParameterList> validParamList = rcp(new ParameterList());
66 
67  validParamList->set< RCP<const FactoryBase> >("Coordinates", Teuchos::null, "Factory for coordinates generation");
68  validParamList->set< RCP<const FactoryBase> >("Aggregates", Teuchos::null, "Factory for coordinates generation");
69  validParamList->set< RCP<const FactoryBase> >("CoarseMap", Teuchos::null, "Generating factory of the coarse map");
70  validParamList->set< int > ("write start", -1, "First level at which coordinates should be written to file");
71  validParamList->set< int > ("write end", -1, "Last level at which coordinates should be written to file");
72  validParamList->set< bool > ("structured aggregation", false, "Flag specifying that the geometric data is transferred for StructuredAggregationFactory");
73  validParamList->set<RCP<const FactoryBase> > ("lCoarseNodesPerDim", Teuchos::null, "Factory providing the local number of nodes per spatial dimensions of the mesh");
74  validParamList->set<RCP<const FactoryBase> > ("numDimensions", Teuchos::null, "Factory providing the number of spatial dimensions of the mesh");
75 
76  return validParamList;
77  }
78 
79  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class DeviceType>
81  static bool isAvailableCoords = false;
82 
83  const ParameterList& pL = GetParameterList();
84  if(pL.get<bool>("structured aggregation") == true) {
85  Input(fineLevel, "lCoarseNodesPerDim");
86  Input(fineLevel, "numDimensions");
87  } else {
88  if (coarseLevel.GetRequestMode() == Level::REQUEST)
89  isAvailableCoords = coarseLevel.IsAvailable("Coordinates", this);
90 
91  if (isAvailableCoords == false) {
92  Input(fineLevel, "Coordinates");
93  Input(fineLevel, "Aggregates");
94  Input(fineLevel, "CoarseMap");
95  }
96  }
97  }
98 
99  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class DeviceType>
101  FactoryMonitor m(*this, "Build", coarseLevel);
102 
104 
105  GetOStream(Runtime0) << "Transferring coordinates" << std::endl;
106 
107  if (coarseLevel.IsAvailable("Coordinates", this)) {
108  GetOStream(Runtime0) << "Reusing coordinates" << std::endl;
109  return;
110  }
111 
112 
113  const ParameterList& pL = GetParameterList();
114  if(pL.get<bool>("structured aggregation") == true) {
115  Array<LO> lCoarseNodesPerDir = Get<Array<LO> >(fineLevel, "lCoarseNodesPerDim");
116  Set<Array<LO> >(coarseLevel, "lNodesPerDim", lCoarseNodesPerDir);
117  int numDimensions = Get<int>(fineLevel, "numDimensions");
118  Set<int>(coarseLevel, "numDimensions", numDimensions);
119 
120  return;
121  }
122 
123  auto aggregates = Get<RCP<Aggregates_kokkos> >(fineLevel, "Aggregates");
124  auto fineCoords = Get<RCP<doubleMultiVector> >(fineLevel, "Coordinates");
125  auto coarseMap = Get<RCP<const Map> > (fineLevel, "CoarseMap");
126 
127  // coarseMap is being used to set up the domain map of tentative P, and therefore, the row map of Ac
128  // Therefore, if we amalgamate coarseMap, logical nodes in the coordinates vector would correspond to
129  // logical blocks in the matrix
130 
131  ArrayView<const GO> elementAList = coarseMap->getLocalElementList();
132  GO indexBase = coarseMap->getIndexBase();
133 
134  LO blkSize = 1;
135  if (rcp_dynamic_cast<const StridedMap>(coarseMap) != Teuchos::null)
136  blkSize = rcp_dynamic_cast<const StridedMap>(coarseMap)->getFixedBlockSize();
137 
138  Array<GO> elementList;
139  ArrayView<const GO> elementListView;
140  if (blkSize == 1) {
141  // Scalar system
142  // No amalgamation required
143  elementListView = elementAList;
144 
145  } else {
146  auto numElements = elementAList.size() / blkSize;
147 
148  elementList.resize(numElements);
149 
150  // Amalgamate the map
151  for (LO i = 0; i < Teuchos::as<LO>(numElements); i++)
152  elementList[i] = (elementAList[i*blkSize]-indexBase)/blkSize + indexBase;
153 
154  elementListView = elementList;
155  }
156 
157  auto uniqueMap = fineCoords->getMap();
158  auto coarseCoordMap = MapFactory::Build(coarseMap->lib(), Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
159  elementListView, indexBase, coarseMap->getComm());
160  RCP<doubleMultiVector> coarseCoords = Xpetra::MultiVectorFactory<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO>::Build(coarseCoordMap, fineCoords->getNumVectors());
161 
162  // Create overlapped fine coordinates to reduce global communication
163  RCP<doubleMultiVector> ghostedCoords = fineCoords;
164  if (aggregates->AggregatesCrossProcessors()) {
165  auto nonUniqueMap = aggregates->GetMap();
166  auto importer = ImportFactory::Build(uniqueMap, nonUniqueMap);
167 
168  ghostedCoords = Xpetra::MultiVectorFactory<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO>::Build(nonUniqueMap, fineCoords->getNumVectors());
169  ghostedCoords->doImport(*fineCoords, *importer, Xpetra::INSERT);
170  }
171 
172  // The good new is that his graph has already been constructed for the
173  // TentativePFactory and was cached in Aggregates. So this is a no-op.
174  auto aggGraph = aggregates->GetGraph();
175  auto numAggs = aggGraph.numRows();
176 
177  auto fineCoordsView = fineCoords ->getDeviceLocalView(Xpetra::Access::ReadOnly);
178  auto coarseCoordsView = coarseCoords->getDeviceLocalView(Xpetra::Access::OverwriteAll);
179 
180  // Fill in coarse coordinates
181  {
182  SubFactoryMonitor m2(*this, "AverageCoords", coarseLevel);
183 
184  const auto dim = fineCoords->getNumVectors();
185 
186  typename AppendTrait<decltype(fineCoordsView), Kokkos::RandomAccess>::type fineCoordsRandomView = fineCoordsView;
187  for (size_t j = 0; j < dim; j++) {
188  Kokkos::parallel_for("MueLu:CoordinatesTransferF:Build:coord", Kokkos::RangePolicy<local_ordinal_type, execution_space>(0, numAggs),
189  KOKKOS_LAMBDA(const LO i) {
190  // A row in this graph represents all node ids in the aggregate
191  // Therefore, averaging is very easy
192 
193  auto aggregate = aggGraph.rowConst(i);
194 
195  typename Teuchos::ScalarTraits<Scalar>::magnitudeType sum = 0.0; // do not use Scalar here (Stokhos)
196  for (size_t colID = 0; colID < static_cast<size_t>(aggregate.length); colID++)
197  sum += fineCoordsRandomView(aggregate(colID),j);
198 
199  coarseCoordsView(i,j) = sum / aggregate.length;
200  });
201  }
202  }
203 
204  Set<RCP<doubleMultiVector> >(coarseLevel, "Coordinates", coarseCoords);
205 
206  int writeStart = pL.get<int>("write start"), writeEnd = pL.get<int>("write end");
207  if (writeStart == 0 && fineLevel.GetLevelID() == 0 && writeStart <= writeEnd) {
208  std::string fileName = "coordinates_before_rebalance_level_" + toString(fineLevel.GetLevelID()) + ".m";
209  Xpetra::IO<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO>::Write(fileName, *fineCoords);
210  }
211  if (writeStart <= coarseLevel.GetLevelID() && coarseLevel.GetLevelID() <= writeEnd) {
212  std::string fileName = "coordinates_before_rebalance_level_" + toString(coarseLevel.GetLevelID()) + ".m";
213  Xpetra::IO<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO>::Write(fileName,*coarseCoords);
214  }
215  }
216 
217 } // namespace MueLu
218 
219 #endif // MUELU_COORDINATESTRANSFER_FACTORY_KOKKOS_DEF_HPP
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
GlobalOrdinal GO
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Timer to be used in factories. Similar to Monitor but with additional timers.
size_type size() const
LocalOrdinal LO
One-liner description of what is happening.
RequestMode GetRequestMode() const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
void resize(size_type new_size, const value_type &x=value_type())
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need&#39;s value has been saved.