MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_IndexManager_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 // Ray Tuminaro (rstumin@sandia.gov)
41 // Luc Berger-Vergiat (lberge@sandia.gov)
42 //
43 // ***********************************************************************
44 //
45 // @HEADER
46 #ifndef MUELU_INDEXMANAGER_DEF_KOKKOS_HPP
47 #define MUELU_INDEXMANAGER_DEF_KOKKOS_HPP
48 
49 #ifdef HAVE_MUELU_KOKKOS_REFACTOR
50 
51 #include <utility>
52 
54 
55 #include <Xpetra_MapFactory.hpp>
56 
57 #include "MueLu_ConfigDefs.hpp"
58 #include "MueLu_Types.hpp"
59 #include "MueLu_BaseClass.hpp"
60 #include "MueLu_Exceptions.hpp"
62 
63 /*****************************************************************************
64 
65 ****************************************************************************/
66 
67 namespace MueLu {
68 
69  template<class LocalOrdinal, class GlobalOrdinal, class Node>
70  IndexManager_kokkos<LocalOrdinal, GlobalOrdinal, Node>::
71  IndexManager_kokkos(const int NumDimensions,
72  const int interpolationOrder,
73  const int MyRank,
74  const ArrayView<const LO> LFineNodesPerDir,
75  const ArrayView<const int> CoarseRate) :
76  myRank(MyRank), coarseRate("coarsening rate"), endRate("endRate"),
77  lFineNodesPerDir("lFineNodesPerDir"), coarseNodesPerDir("lFineNodesPerDir") {
78 
79  RCP<Teuchos::FancyOStream> out;
80  if(const char* dbg = std::getenv("MUELU_INDEXMANAGER_DEBUG")) {
81  out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
82  out->setShowAllFrontMatter(false).setShowProcRank(true);
83  } else {
84  out = Teuchos::getFancyOStream(rcp(new Teuchos::oblackholestream()));
85  }
86 
87  setupIM(NumDimensions, interpolationOrder, CoarseRate, LFineNodesPerDir);
88 
89  *out << "Done setting up the IndexManager" << std::endl;
90 
91  computeMeshParameters();
92 
93  *out << "Computed Mesh Parameters" << std::endl;
94 
95  } // IndexManager_kokkos Constructor
96 
97  template<class LocalOrdinal, class GlobalOrdinal, class Node>
98  void IndexManager_kokkos<LocalOrdinal, GlobalOrdinal, Node>::
99  setupIM(const int NumDimensions, const int interpolationOrder,
100  const ArrayView<const int> CoarseRate, const ArrayView<const LO> LFineNodesPerDir) {
101 
102  numDimensions = NumDimensions;
103  interpolationOrder_ = interpolationOrder;
104 
105  TEUCHOS_TEST_FOR_EXCEPTION((LFineNodesPerDir.size() != 3)
106  && (LFineNodesPerDir.size() != numDimensions),
107  Exceptions::RuntimeError,
108  "LFineNodesPerDir has to be of size 3 or of size numDimensions!");
109 
110  typename Kokkos::View<LO[3], memory_space>::HostMirror lFineNodesPerDir_h = Kokkos::create_mirror_view(lFineNodesPerDir);
111  Kokkos::deep_copy(lFineNodesPerDir_h, lFineNodesPerDir);
112  typename Kokkos::View<LO[3], memory_space>::HostMirror coarseRate_h = Kokkos::create_mirror_view(coarseRate);
113  Kokkos::deep_copy(coarseRate_h, coarseRate);
114 
115  // Load coarse rate, being careful about formating
116  // Also load lFineNodesPerDir
117  for(int dim = 0; dim < 3; ++dim) {
118  if(dim < getNumDimensions()) {
119  lFineNodesPerDir_h(dim) = LFineNodesPerDir[dim];
120  if(CoarseRate.size() == 1) {
121  coarseRate_h(dim) = CoarseRate[0];
122  } else if(CoarseRate.size() == getNumDimensions()) {
123  coarseRate_h(dim) = CoarseRate[dim];
124  }
125  } else {
126  lFineNodesPerDir_h(dim) = 1;
127  coarseRate_h(dim) = 1;
128  }
129  }
130 
131  Kokkos::deep_copy(lFineNodesPerDir, lFineNodesPerDir_h);
132  Kokkos::deep_copy(coarseRate, coarseRate_h);
133 
134  } // setupIM
135 
136  template<class LocalOrdinal, class GlobalOrdinal, class Node>
137  void IndexManager_kokkos<LocalOrdinal, GlobalOrdinal, Node>::computeMeshParameters() {
138 
139  RCP<Teuchos::FancyOStream> out;
140  if(const char* dbg = std::getenv("MUELU_INDEXMANAGER_DEBUG")) {
141  out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
142  out->setShowAllFrontMatter(false).setShowProcRank(true);
143  } else {
144  out = Teuchos::getFancyOStream(rcp(new Teuchos::oblackholestream()));
145  }
146 
147  typename Kokkos::View<int[3], memory_space>::HostMirror coarseRate_h = Kokkos::create_mirror_view(coarseRate);
148  typename Kokkos::View<int[3], memory_space>::HostMirror endRate_h = Kokkos::create_mirror_view(endRate);
149 
150 
151  typename Kokkos::View<LO[3], memory_space>::HostMirror lFineNodesPerDir_h = Kokkos::create_mirror_view(lFineNodesPerDir);
152  typename Kokkos::View<LO[3], memory_space>::HostMirror coarseNodesPerDir_h = Kokkos::create_mirror_view(coarseNodesPerDir);
153  Kokkos::deep_copy(lFineNodesPerDir_h, lFineNodesPerDir);
154 
155  lNumFineNodes10 = lFineNodesPerDir_h(1)*lFineNodesPerDir_h(0);
156  lNumFineNodes = lFineNodesPerDir_h(2)*lNumFineNodes10;
157  for(int dim = 0; dim < 3; ++dim) {
158  if(dim < numDimensions) {
159  endRate_h(dim) = (lFineNodesPerDir_h(dim) - 1) % coarseRate_h(dim);
160  if(endRate_h(dim) == 0) {endRate_h(dim) = coarseRate_h(dim);}
161 
162  } else { // Default value for dim >= numDimensions
163  endRate_h(dim) = 1;
164  }
165  }
166 
167  *out << "lFineNodesPerDir: {" << lFineNodesPerDir_h(0) << ", " << lFineNodesPerDir_h(1) << ", "
168  << lFineNodesPerDir_h(2) << "}" << std::endl;
169  *out << "endRate: {" << endRate_h(0) << ", " << endRate_h(1) << ", "
170  << endRate_h(2) << "}" << std::endl;
171 
172  // Here one element can represent either the degenerate case of one node or the more general
173  // case of two nodes, i.e. x---x is a 1D element with two nodes and x is a 1D element with
174  // one node. This helps generating a 3D space from tensorial products...
175  // A good way to handle this would be to generalize the algorithm to take into account the
176  // discretization order used in each direction, at least in the FEM sense, since a 0 degree
177  // discretization will have a unique node per element. This way 1D discretization can be
178  // viewed as a 3D problem with one 0 degree element in the y direction and one 0 degre
179  // element in the z direction.
180  // !!! Operations below are aftecting both local and global values that have two !!!
181  // different orientations. Orientations can be interchanged using mapDirG2L and mapDirL2G.
182  // coarseRate, endRate and offsets are in the global basis, as well as all the variables
183  // starting with a g.
184  // !!! while the variables starting with an l are in the local basis. !!!
185  for(int dim = 0; dim < 3; ++dim) {
186  if(dim < numDimensions) {
187  // Check whether the partition includes the "end" of the mesh which means that endRate
188  // will apply. Also make sure that endRate is not 0 which means that the mesh does not
189  // require a particular treatment at the boundaries.
190  coarseNodesPerDir_h(dim) = (lFineNodesPerDir_h(dim) - endRate_h(dim) - 1)
191  / coarseRate_h(dim) + 2;
192 
193  } else { // Default value for dim >= numDimensions
194  // endRate[dim] = 1;
195  coarseNodesPerDir_h(dim) = 1;
196  } // if (dim < numDimensions)
197 
198  // This would happen if the rank does not own any nodes but in that case a subcommunicator
199  // should be used so this should really not be a concern.
200  if(lFineNodesPerDir_h(dim) < 1) {coarseNodesPerDir_h(dim) = 0;}
201  } // Loop for dim=0:3
202 
203  // Compute cummulative values
204  numCoarseNodes10 = coarseNodesPerDir_h(0)*coarseNodesPerDir_h(1);
205  numCoarseNodes = numCoarseNodes10*coarseNodesPerDir_h(2);
206 
207  *out << "coarseNodesPerDir: {" << coarseNodesPerDir_h(0) << ", "
208  << coarseNodesPerDir_h(1) << ", " << coarseNodesPerDir_h(2) << "}" << std::endl;
209  *out << "numCoarseNodes=" << numCoarseNodes << std::endl;
210 
211  // Copy Host data to Device.
212  Kokkos::deep_copy(coarseRate, coarseRate_h);
213  Kokkos::deep_copy(endRate, endRate_h);
214  Kokkos::deep_copy(lFineNodesPerDir, lFineNodesPerDir_h);
215  Kokkos::deep_copy(coarseNodesPerDir, coarseNodesPerDir_h);
216  }
217 
218  template<class LocalOrdinal, class GlobalOrdinal, class Node>
219  Array<LocalOrdinal> IndexManager_kokkos<LocalOrdinal, GlobalOrdinal, Node>::
220  getCoarseNodesPerDirArray() const {
221  typename LOTupleView::HostMirror coarseNodesPerDir_h = Kokkos::create_mirror_view(coarseNodesPerDir);
222  Array<LO> coarseNodesPerDirArray(3);
223 
224  for(int dim = 0; dim < 3; ++dim) {
225  coarseNodesPerDirArray[dim] = coarseNodesPerDir_h(dim);
226  }
227 
228  return coarseNodesPerDirArray;
229  } // getCoarseNodesData
230 
231 } //namespace MueLu
232 
233 #endif // HAVE_MUELU_KOKKOS_REFACTOR
234 #define MUELU_INDEXMANAGER_KOKKOS_SHORT
235 #endif // MUELU_INDEXMANAGER_DEF_KOKKOS_HPP
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void deep_copy(const View< DT, DP...> &dst, typename ViewTraits< DT, DP...>::const_value_type &value, typename std::enable_if< std::is_same< typename ViewTraits< DT, DP...>::specialize, void >::value >::type *=0)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)