MueLu Version of the Day
Loading...
Searching...
No Matches
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_CoarseMapFactory_kokkos.hpp"
58#include "MueLu_Level.hpp"
59#include "MueLu_Monitor.hpp"
60#include "MueLu_Utilities_kokkos.hpp"
61
62namespace MueLu {
63
64 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class DeviceType>
66 RCP<ParameterList> validParamList = rcp(new ParameterList());
67
68 validParamList->set< RCP<const FactoryBase> >("Coordinates", Teuchos::null, "Factory for coordinates generation");
69 validParamList->set< RCP<const FactoryBase> >("Aggregates", Teuchos::null, "Factory for coordinates generation");
70 validParamList->set< RCP<const FactoryBase> >("CoarseMap", Teuchos::null, "Generating factory of the coarse map");
71 validParamList->set< int > ("write start", -1, "First level at which coordinates should be written to file");
72 validParamList->set< int > ("write end", -1, "Last level at which coordinates should be written to file");
73 validParamList->set< bool > ("structured aggregation", false, "Flag specifying that the geometric data is transferred for StructuredAggregationFactory");
74 validParamList->set<RCP<const FactoryBase> > ("lCoarseNodesPerDim", Teuchos::null, "Factory providing the local number of nodes per spatial dimensions of the mesh");
75 validParamList->set<RCP<const FactoryBase> > ("numDimensions", Teuchos::null, "Factory providing the number of spatial dimensions of the mesh");
76
77 return validParamList;
78 }
79
80 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class DeviceType>
82 static bool isAvailableCoords = false;
83
84 const ParameterList& pL = GetParameterList();
85 if(pL.get<bool>("structured aggregation") == true) {
86 Input(fineLevel, "lCoarseNodesPerDim");
87 Input(fineLevel, "numDimensions");
88 } else {
89 if (coarseLevel.GetRequestMode() == Level::REQUEST)
90 isAvailableCoords = coarseLevel.IsAvailable("Coordinates", this);
91
92 if (isAvailableCoords == false) {
93 Input(fineLevel, "Coordinates");
94 Input(fineLevel, "Aggregates");
95 Input(fineLevel, "CoarseMap");
96 }
97 }
98 }
99
100 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class DeviceType>
102 FactoryMonitor m(*this, "Build", coarseLevel);
103
104 typedef Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO> doubleMultiVector;
105
106 GetOStream(Runtime0) << "Transferring coordinates" << std::endl;
107
108 if (coarseLevel.IsAvailable("Coordinates", this)) {
109 GetOStream(Runtime0) << "Reusing coordinates" << std::endl;
110 return;
111 }
112
113
114 const ParameterList& pL = GetParameterList();
115 if(pL.get<bool>("structured aggregation") == true) {
116 Array<LO> lCoarseNodesPerDir = Get<Array<LO> >(fineLevel, "lCoarseNodesPerDim");
117 Set<Array<LO> >(coarseLevel, "lNodesPerDim", lCoarseNodesPerDir);
118 int numDimensions = Get<int>(fineLevel, "numDimensions");
119 Set<int>(coarseLevel, "numDimensions", numDimensions);
120
121 return;
122 }
123
124 auto aggregates = Get<RCP<Aggregates_kokkos> >(fineLevel, "Aggregates");
125 auto fineCoords = Get<RCP<doubleMultiVector> >(fineLevel, "Coordinates");
126 auto coarseMap = Get<RCP<const Map> > (fineLevel, "CoarseMap");
127
128 // coarseMap is being used to set up the domain map of tentative P, and therefore, the row map of Ac
129 // Therefore, if we amalgamate coarseMap, logical nodes in the coordinates vector would correspond to
130 // logical blocks in the matrix
131
132 ArrayView<const GO> elementAList = coarseMap->getLocalElementList();
133 GO indexBase = coarseMap->getIndexBase();
134
135 LO blkSize = 1;
136 if (rcp_dynamic_cast<const StridedMap>(coarseMap) != Teuchos::null)
137 blkSize = rcp_dynamic_cast<const StridedMap>(coarseMap)->getFixedBlockSize();
138
139 Array<GO> elementList;
140 ArrayView<const GO> elementListView;
141 if (blkSize == 1) {
142 // Scalar system
143 // No amalgamation required
144 elementListView = elementAList;
145
146 } else {
147 auto numElements = elementAList.size() / blkSize;
148
149 elementList.resize(numElements);
150
151 // Amalgamate the map
152 for (LO i = 0; i < Teuchos::as<LO>(numElements); i++)
153 elementList[i] = (elementAList[i*blkSize]-indexBase)/blkSize + indexBase;
154
155 elementListView = elementList;
156 }
157
158 auto uniqueMap = fineCoords->getMap();
159 auto coarseCoordMap = MapFactory::Build(coarseMap->lib(), Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
160 elementListView, indexBase, coarseMap->getComm());
161 RCP<doubleMultiVector> coarseCoords = Xpetra::MultiVectorFactory<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO>::Build(coarseCoordMap, fineCoords->getNumVectors());
162
163 // Create overlapped fine coordinates to reduce global communication
164 RCP<doubleMultiVector> ghostedCoords = fineCoords;
165 if (aggregates->AggregatesCrossProcessors()) {
166 auto nonUniqueMap = aggregates->GetMap();
167 auto importer = ImportFactory::Build(uniqueMap, nonUniqueMap);
168
169 ghostedCoords = Xpetra::MultiVectorFactory<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO>::Build(nonUniqueMap, fineCoords->getNumVectors());
170 ghostedCoords->doImport(*fineCoords, *importer, Xpetra::INSERT);
171 }
172
173 // The good new is that his graph has already been constructed for the
174 // TentativePFactory and was cached in Aggregates. So this is a no-op.
175 auto aggGraph = aggregates->GetGraph();
176 auto numAggs = aggGraph.numRows();
177
178 auto fineCoordsView = fineCoords ->getDeviceLocalView(Xpetra::Access::ReadOnly);
179 auto coarseCoordsView = coarseCoords->getDeviceLocalView(Xpetra::Access::OverwriteAll);
180
181 // Fill in coarse coordinates
182 {
183 SubFactoryMonitor m2(*this, "AverageCoords", coarseLevel);
184
185 const auto dim = fineCoords->getNumVectors();
186
187 typename AppendTrait<decltype(fineCoordsView), Kokkos::RandomAccess>::type fineCoordsRandomView = fineCoordsView;
188 for (size_t j = 0; j < dim; j++) {
189 Kokkos::parallel_for("MueLu:CoordinatesTransferF:Build:coord", Kokkos::RangePolicy<local_ordinal_type, execution_space>(0, numAggs),
190 KOKKOS_LAMBDA(const LO i) {
191 // A row in this graph represents all node ids in the aggregate
192 // Therefore, averaging is very easy
193
194 auto aggregate = aggGraph.rowConst(i);
195
196 typename Teuchos::ScalarTraits<Scalar>::magnitudeType sum = 0.0; // do not use Scalar here (Stokhos)
197 for (size_t colID = 0; colID < static_cast<size_t>(aggregate.length); colID++)
198 sum += fineCoordsRandomView(aggregate(colID),j);
199
200 coarseCoordsView(i,j) = sum / aggregate.length;
201 });
202 }
203 }
204
205 Set<RCP<doubleMultiVector> >(coarseLevel, "Coordinates", coarseCoords);
206
207 int writeStart = pL.get<int>("write start"), writeEnd = pL.get<int>("write end");
208 if (writeStart == 0 && fineLevel.GetLevelID() == 0 && writeStart <= writeEnd) {
209 std::string fileName = "coordinates_before_rebalance_level_" + toString(fineLevel.GetLevelID()) + ".m";
210 Xpetra::IO<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO>::Write(fileName, *fineCoords);
211 }
212 if (writeStart <= coarseLevel.GetLevelID() && coarseLevel.GetLevelID() <= writeEnd) {
213 std::string fileName = "coordinates_before_rebalance_level_" + toString(coarseLevel.GetLevelID()) + ".m";
214 Xpetra::IO<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO>::Write(fileName,*coarseCoords);
215 }
216 }
217
218} // namespace MueLu
219
220#endif // MUELU_COORDINATESTRANSFER_FACTORY_KOKKOS_DEF_HPP
Timer to be used in factories. Similar to Monitor but with additional timers.
Class that holds all level-specific information.
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
int GetLevelID() const
Return level number.
RequestMode GetRequestMode() const
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
Namespace for MueLu classes and methods.
@ Runtime0
One-liner description of what is happening.
std::string toString(const T &what)
Little helper function to convert non-string types to strings.