MueLu Version of the Day
Loading...
Searching...
No Matches
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#include <utility>
50
51#include "Teuchos_OrdinalTraits.hpp"
52
53#include <Xpetra_MapFactory.hpp>
54
55#include "MueLu_ConfigDefs.hpp"
56#include "MueLu_Types.hpp"
57#include "MueLu_BaseClass.hpp"
58#include "MueLu_Exceptions.hpp"
60
61/*****************************************************************************
62
63****************************************************************************/
64
65namespace MueLu {
66
67 template<class LocalOrdinal, class GlobalOrdinal, class Node>
69 IndexManager_kokkos(const int NumDimensions,
70 const int interpolationOrder,
71 const int MyRank,
72 const ArrayView<const LO> LFineNodesPerDir,
73 const ArrayView<const int> CoarseRate) :
74 myRank(MyRank), coarseRate("coarsening rate"), endRate("endRate"),
75 lFineNodesPerDir("lFineNodesPerDir"), coarseNodesPerDir("lFineNodesPerDir") {
76
77 RCP<Teuchos::FancyOStream> out;
78 if(const char* dbg = std::getenv("MUELU_INDEXMANAGER_DEBUG")) {
79 out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
80 out->setShowAllFrontMatter(false).setShowProcRank(true);
81 } else {
82 out = Teuchos::getFancyOStream(rcp(new Teuchos::oblackholestream()));
83 }
84
85 setupIM(NumDimensions, interpolationOrder, CoarseRate, LFineNodesPerDir);
86
87 *out << "Done setting up the IndexManager" << std::endl;
88
90
91 *out << "Computed Mesh Parameters" << std::endl;
92
93 } // IndexManager_kokkos Constructor
94
95 template<class LocalOrdinal, class GlobalOrdinal, class Node>
97 setupIM(const int NumDimensions, const int interpolationOrder,
98 const ArrayView<const int> CoarseRate, const ArrayView<const LO> LFineNodesPerDir) {
99
100 numDimensions = NumDimensions;
101 interpolationOrder_ = interpolationOrder;
102
103 TEUCHOS_TEST_FOR_EXCEPTION((LFineNodesPerDir.size() != 3)
104 && (LFineNodesPerDir.size() != numDimensions),
106 "LFineNodesPerDir has to be of size 3 or of size numDimensions!");
107
108 typename Kokkos::View<LO[3], device_type>::HostMirror lFineNodesPerDir_h = Kokkos::create_mirror_view(lFineNodesPerDir);
109 Kokkos::deep_copy(lFineNodesPerDir_h, lFineNodesPerDir);
110 typename Kokkos::View<LO[3], device_type>::HostMirror coarseRate_h = Kokkos::create_mirror_view(coarseRate);
111 Kokkos::deep_copy(coarseRate_h, coarseRate);
112
113 // Load coarse rate, being careful about formating
114 // Also load lFineNodesPerDir
115 for(int dim = 0; dim < 3; ++dim) {
116 if(dim < getNumDimensions()) {
117 lFineNodesPerDir_h(dim) = LFineNodesPerDir[dim];
118 if(CoarseRate.size() == 1) {
119 coarseRate_h(dim) = CoarseRate[0];
120 } else if(CoarseRate.size() == getNumDimensions()) {
121 coarseRate_h(dim) = CoarseRate[dim];
122 }
123 } else {
124 lFineNodesPerDir_h(dim) = 1;
125 coarseRate_h(dim) = 1;
126 }
127 }
128
129 Kokkos::deep_copy(lFineNodesPerDir, lFineNodesPerDir_h);
130 Kokkos::deep_copy(coarseRate, coarseRate_h);
131
132 } // setupIM
133
134 template<class LocalOrdinal, class GlobalOrdinal, class Node>
136
137 RCP<Teuchos::FancyOStream> out;
138 if(const char* dbg = std::getenv("MUELU_INDEXMANAGER_DEBUG")) {
139 out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
140 out->setShowAllFrontMatter(false).setShowProcRank(true);
141 } else {
142 out = Teuchos::getFancyOStream(rcp(new Teuchos::oblackholestream()));
143 }
144
145 typename Kokkos::View<int[3], device_type>::HostMirror coarseRate_h = Kokkos::create_mirror_view(coarseRate);
146 typename Kokkos::View<int[3], device_type>::HostMirror endRate_h = Kokkos::create_mirror_view(endRate);
147
148
149 typename Kokkos::View<LO[3], device_type>::HostMirror lFineNodesPerDir_h = Kokkos::create_mirror_view(lFineNodesPerDir);
150 typename Kokkos::View<LO[3], device_type>::HostMirror coarseNodesPerDir_h = Kokkos::create_mirror_view(coarseNodesPerDir);
151 Kokkos::deep_copy(lFineNodesPerDir_h, lFineNodesPerDir);
152 Kokkos::deep_copy(coarseRate_h, coarseRate);
153
154 lNumFineNodes10 = lFineNodesPerDir_h(1)*lFineNodesPerDir_h(0);
155 lNumFineNodes = lFineNodesPerDir_h(2)*lNumFineNodes10;
156 for(int dim = 0; dim < 3; ++dim) {
157 if(dim < numDimensions) {
158 endRate_h(dim) = (lFineNodesPerDir_h(dim) - 1) % coarseRate_h(dim);
159 if(endRate_h(dim) == 0) {endRate_h(dim) = coarseRate_h(dim);}
160
161 } else { // Default value for dim >= numDimensions
162 endRate_h(dim) = 1;
163 }
164 }
165
166 *out << "lFineNodesPerDir: {" << lFineNodesPerDir_h(0) << ", " << lFineNodesPerDir_h(1) << ", "
167 << lFineNodesPerDir_h(2) << "}" << std::endl;
168 *out << "endRate: {" << endRate_h(0) << ", " << endRate_h(1) << ", "
169 << endRate_h(2) << "}" << std::endl;
170
171 // Here one element can represent either the degenerate case of one node or the more general
172 // case of two nodes, i.e. x---x is a 1D element with two nodes and x is a 1D element with
173 // one node. This helps generating a 3D space from tensorial products...
174 // A good way to handle this would be to generalize the algorithm to take into account the
175 // discretization order used in each direction, at least in the FEM sense, since a 0 degree
176 // discretization will have a unique node per element. This way 1D discretization can be
177 // viewed as a 3D problem with one 0 degree element in the y direction and one 0 degre
178 // element in the z direction.
179 // !!! Operations below are aftecting both local and global values that have two !!!
180 // different orientations. Orientations can be interchanged using mapDirG2L and mapDirL2G.
181 // coarseRate, endRate and offsets are in the global basis, as well as all the variables
182 // starting with a g.
183 // !!! while the variables starting with an l are in the local basis. !!!
184 for(int dim = 0; dim < 3; ++dim) {
185 if(dim < numDimensions) {
186 // Check whether the partition includes the "end" of the mesh which means that endRate
187 // will apply. Also make sure that endRate is not 0 which means that the mesh does not
188 // require a particular treatment at the boundaries.
189 coarseNodesPerDir_h(dim) = (lFineNodesPerDir_h(dim) - endRate_h(dim) - 1)
190 / coarseRate_h(dim) + 2;
191
192 } else { // Default value for dim >= numDimensions
193 // endRate[dim] = 1;
194 coarseNodesPerDir_h(dim) = 1;
195 } // if (dim < numDimensions)
196
197 // This would happen if the rank does not own any nodes but in that case a subcommunicator
198 // should be used so this should really not be a concern.
199 if(lFineNodesPerDir_h(dim) < 1) {coarseNodesPerDir_h(dim) = 0;}
200 } // Loop for dim=0:3
201
202 // Compute cummulative values
203 numCoarseNodes10 = coarseNodesPerDir_h(0)*coarseNodesPerDir_h(1);
204 numCoarseNodes = numCoarseNodes10*coarseNodesPerDir_h(2);
205
206 *out << "coarseNodesPerDir: {" << coarseNodesPerDir_h(0) << ", "
207 << coarseNodesPerDir_h(1) << ", " << coarseNodesPerDir_h(2) << "}" << std::endl;
208 *out << "numCoarseNodes=" << numCoarseNodes << std::endl;
209
210 // Copy Host data to Device.
211 Kokkos::deep_copy(coarseRate, coarseRate_h);
212 Kokkos::deep_copy(endRate, endRate_h);
213 Kokkos::deep_copy(lFineNodesPerDir, lFineNodesPerDir_h);
214 Kokkos::deep_copy(coarseNodesPerDir, coarseNodesPerDir_h);
215 }
216
217 template<class LocalOrdinal, class GlobalOrdinal, class Node>
220 typename LOTupleView::HostMirror coarseNodesPerDir_h = Kokkos::create_mirror_view(coarseNodesPerDir);
221 Kokkos::deep_copy(coarseNodesPerDir_h, 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#define MUELU_INDEXMANAGER_KOKKOS_SHORT
234#endif // MUELU_INDEXMANAGER_DEF_KOKKOS_HPP
Exception throws to report errors in the internal logical of the program.
IndexManager_kokkos()=default
Default constructor, return empty object.
void setupIM(const int NumDimensions, const int interpolationOrder, const ArrayView< const int > coarseRate, const ArrayView< const LO > LFineNodesPerDir)
Common setup pattern used for all the different types of undelying mesh.
Namespace for MueLu classes and methods.