MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_CoordinatesTransferFactory_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_DEF_HPP
47#define MUELU_COORDINATESTRANSFER_FACTORY_DEF_HPP
48
49#include "Xpetra_ImportFactory.hpp"
50#include "Xpetra_MultiVectorFactory.hpp"
51#include "Xpetra_MapFactory.hpp"
52#include "Xpetra_IO.hpp"
53
54#include "MueLu_CoarseMapFactory.hpp"
55#include "MueLu_Aggregates.hpp"
57//#include "MueLu_Utilities.hpp"
58
59#include "MueLu_Level.hpp"
60#include "MueLu_Monitor.hpp"
61
62namespace MueLu {
63
64 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
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<bool> ("structured aggregation", false, "Flag specifying that the geometric data is transferred for StructuredAggregationFactory");
72 validParamList->set<bool> ("aggregation coupled", false, "Flag specifying if the aggregation algorithm was used in coupled mode.");
73 validParamList->set<bool> ("Geometric", false, "Flag specifying that the coordinates are transferred for GeneralGeometricPFactory");
74 validParamList->set<RCP<const FactoryBase> >("coarseCoordinates", Teuchos::null, "Factory for coarse coordinates generation");
75 validParamList->set<RCP<const FactoryBase> >("gCoarseNodesPerDim", Teuchos::null, "Factory providing the global number of nodes per spatial dimensions of the mesh");
76 validParamList->set<RCP<const FactoryBase> >("lCoarseNodesPerDim", Teuchos::null, "Factory providing the local number of nodes per spatial dimensions of the mesh");
77 validParamList->set<RCP<const FactoryBase> >("numDimensions" , Teuchos::null, "Factory providing the number of spatial dimensions of the mesh");
78 validParamList->set<int> ("write start", -1, "first level at which coordinates should be written to file");
79 validParamList->set<int> ("write end", -1, "last level at which coordinates should be written to file");
80 validParamList->set<bool> ("hybrid aggregation", false, "Flag specifying that hybrid aggregation data is transfered for HybridAggregationFactory");
81 validParamList->set<RCP<const FactoryBase> >("aggregationRegionTypeCoarse", Teuchos::null, "Factory indicating what aggregation type is to be used on the coarse level of the region");
82 validParamList->set<bool> ("interface aggregation", false, "Flag specifying that interface aggregation data is transfered for HybridAggregationFactory");
83 validParamList->set<RCP<const FactoryBase> >("coarseInterfacesDimensions", Teuchos::null, "Factory providing coarseInterfacesDimensions");
84 validParamList->set<RCP<const FactoryBase> >("nodeOnCoarseInterface", Teuchos::null, "Factory providing nodeOnCoarseInterface");
85
86
87 return validParamList;
88 }
89
90 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
92 static bool isAvailableCoords = false;
93
94 const ParameterList& pL = GetParameterList();
95 if(pL.get<bool>("structured aggregation") == true) {
96 if(pL.get<bool>("aggregation coupled") == true) {
97 Input(fineLevel, "gCoarseNodesPerDim");
98 }
99 Input(fineLevel, "lCoarseNodesPerDim");
100 Input(fineLevel, "numDimensions");
101 } else if(pL.get<bool>("Geometric") == true) {
102 Input(coarseLevel, "coarseCoordinates");
103 Input(coarseLevel, "gCoarseNodesPerDim");
104 Input(coarseLevel, "lCoarseNodesPerDim");
105 } else if(pL.get<bool>("hybrid aggregation") == true) {
106 Input(fineLevel, "aggregationRegionTypeCoarse");
107 Input(fineLevel, "lCoarseNodesPerDim");
108 Input(fineLevel, "numDimensions");
109 if(pL.get<bool>("interface aggregation") == true) {
110 Input(fineLevel, "coarseInterfacesDimensions");
111 Input(fineLevel, "nodeOnCoarseInterface");
112 }
113 } else {
114 if (coarseLevel.GetRequestMode() == Level::REQUEST)
115 isAvailableCoords = coarseLevel.IsAvailable("Coordinates", this);
116
117 if (isAvailableCoords == false) {
118 Input(fineLevel, "Coordinates");
119 Input(fineLevel, "Aggregates");
120 Input(fineLevel, "CoarseMap");
121 }
122 }
123 }
124
125 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
127 FactoryMonitor m(*this, "Build", coarseLevel);
128
129 using xdMV = Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType,LO,GO,NO>;
130
131 GetOStream(Runtime0) << "Transferring coordinates" << std::endl;
132
133 int numDimensions;
134 RCP<xdMV> coarseCoords;
135 RCP<xdMV> fineCoords;
136 Array<GO> gCoarseNodesPerDir;
137 Array<LO> lCoarseNodesPerDir;
138
139 const ParameterList& pL = GetParameterList();
140
141 if(pL.get<bool>("hybrid aggregation") == true) {
142 std::string regionType = Get<std::string>(fineLevel,"aggregationRegionTypeCoarse");
143 numDimensions = Get<int>(fineLevel, "numDimensions");
144 lCoarseNodesPerDir = Get<Array<LO> >(fineLevel, "lCoarseNodesPerDim");
145 Set<std::string>(coarseLevel, "aggregationRegionType", regionType);
146 Set<int> (coarseLevel, "numDimensions", numDimensions);
147 Set<Array<LO> > (coarseLevel, "lNodesPerDim", lCoarseNodesPerDir);
148
149 if((pL.get<bool>("interface aggregation") == true) && (regionType == "uncoupled")) {
150 Array<LO> coarseInterfacesDimensions = Get<Array<LO> >(fineLevel, "coarseInterfacesDimensions");
151 Array<LO> nodeOnCoarseInterface = Get<Array<LO> >(fineLevel, "nodeOnCoarseInterface");
152 Set<Array<LO> >(coarseLevel, "interfacesDimensions", coarseInterfacesDimensions);
153 Set<Array<LO> >(coarseLevel, "nodeOnInterface", nodeOnCoarseInterface);
154 }
155
156 } else if(pL.get<bool>("structured aggregation") == true) {
157 if(pL.get<bool>("aggregation coupled") == true) {
158 gCoarseNodesPerDir = Get<Array<GO> >(fineLevel, "gCoarseNodesPerDim");
159 Set<Array<GO> >(coarseLevel, "gNodesPerDim", gCoarseNodesPerDir);
160 }
161 lCoarseNodesPerDir = Get<Array<LO> >(fineLevel, "lCoarseNodesPerDim");
162 Set<Array<LO> >(coarseLevel, "lNodesPerDim", lCoarseNodesPerDir);
163 numDimensions = Get<int>(fineLevel, "numDimensions");
164 Set<int>(coarseLevel, "numDimensions", numDimensions);
165
166 } else if(pL.get<bool>("Geometric") == true) {
167 coarseCoords = Get<RCP<xdMV> >(coarseLevel, "coarseCoordinates");
168 gCoarseNodesPerDir = Get<Array<GO> >(coarseLevel, "gCoarseNodesPerDim");
169 lCoarseNodesPerDir = Get<Array<LO> >(coarseLevel, "lCoarseNodesPerDim");
170 Set<Array<GO> >(coarseLevel, "gNodesPerDim", gCoarseNodesPerDir);
171 Set<Array<LO> >(coarseLevel, "lNodesPerDim", lCoarseNodesPerDir);
172
173 Set<RCP<xdMV> >(coarseLevel, "Coordinates", coarseCoords);
174
175 } else {
176 if (coarseLevel.IsAvailable("Coordinates", this)) {
177 GetOStream(Runtime0) << "Reusing coordinates" << std::endl;
178 return;
179 }
180
181 RCP<Aggregates> aggregates = Get< RCP<Aggregates> > (fineLevel, "Aggregates");
182 fineCoords = Get< RCP<xdMV> >(fineLevel, "Coordinates");
183 RCP<const Map> coarseMap = Get< RCP<const Map> > (fineLevel, "CoarseMap");
184
185 // coarseMap is being used to set up the domain map of tentative P, and therefore, the row map of Ac
186 // Therefore, if we amalgamate coarseMap, logical nodes in the coordinates vector would correspond to
187 // logical blocks in the matrix
188
189 ArrayView<const GO> elementAList = coarseMap->getLocalElementList();
190
191 LO blkSize = 1;
192 if (rcp_dynamic_cast<const StridedMap>(coarseMap) != Teuchos::null)
193 blkSize = rcp_dynamic_cast<const StridedMap>(coarseMap)->getFixedBlockSize();
194
195 GO indexBase = coarseMap->getIndexBase();
196 size_t numElements = elementAList.size() / blkSize;
197 Array<GO> elementList(numElements);
198
199 // Amalgamate the map
200 for (LO i = 0; i < Teuchos::as<LO>(numElements); i++)
201 elementList[i] = (elementAList[i*blkSize]-indexBase)/blkSize + indexBase;
202
203 RCP<const Map> uniqueMap = fineCoords->getMap();
204 RCP<const Map> coarseCoordMap = MapFactory ::Build(coarseMap->lib(), Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(), elementList, indexBase, coarseMap->getComm());
205 coarseCoords = Xpetra::MultiVectorFactory<typename Teuchos::ScalarTraits<Scalar>::coordinateType,LO,GO,NO>::Build(coarseCoordMap, fineCoords->getNumVectors());
206
207 // Create overlapped fine coordinates to reduce global communication
208 RCP<xdMV> ghostedCoords = fineCoords;
209 if (aggregates->AggregatesCrossProcessors()) {
210 RCP<const Map> nonUniqueMap = aggregates->GetMap();
211 RCP<const Import> importer = ImportFactory::Build(uniqueMap, nonUniqueMap);
212
213 ghostedCoords = Xpetra::MultiVectorFactory<typename Teuchos::ScalarTraits<Scalar>::coordinateType,LO,GO,NO>::Build(nonUniqueMap, fineCoords->getNumVectors());
214 ghostedCoords->doImport(*fineCoords, *importer, Xpetra::INSERT);
215 }
216
217 // Get some info about aggregates
218 int myPID = uniqueMap->getComm()->getRank();
219 LO numAggs = aggregates->GetNumAggregates();
220 ArrayRCP<LO> aggSizes = aggregates->ComputeAggregateSizes();
221 const ArrayRCP<const LO> vertex2AggID = aggregates->GetVertex2AggId()->getData(0);
222 const ArrayRCP<const LO> procWinner = aggregates->GetProcWinner()->getData(0);
223
224 // Fill in coarse coordinates
225 for (size_t j = 0; j < fineCoords->getNumVectors(); j++) {
226 ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> fineCoordsData = ghostedCoords->getData(j);
227 ArrayRCP<typename Teuchos::ScalarTraits<Scalar>::coordinateType> coarseCoordsData = coarseCoords->getDataNonConst(j);
228
229 for (LO lnode = 0; lnode < vertex2AggID.size(); lnode++) {
230 if (procWinner[lnode] == myPID &&
231 lnode < vertex2AggID.size() &&
232 lnode < fineCoordsData.size() && // TAW do not access off-processor coordinates
233 vertex2AggID[lnode] < coarseCoordsData.size() &&
234 Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<Scalar>::coordinateType>::isnaninf(fineCoordsData[lnode]) == false) {
235 coarseCoordsData[vertex2AggID[lnode]] += fineCoordsData[lnode];
236 }
237 }
238 for (LO agg = 0; agg < numAggs; agg++) {
239 coarseCoordsData[agg] /= aggSizes[agg];
240 }
241 }
242
243 Set<RCP<xdMV> >(coarseLevel, "Coordinates", coarseCoords);
244 } // if pL.get<bool>("Geometric") == true
245
246 int writeStart = pL.get<int>("write start"), writeEnd = pL.get<int>("write end");
247 if (writeStart == 0 && fineLevel.GetLevelID() == 0 && writeStart <= writeEnd) {
248 std::ostringstream buf;
249 buf << fineLevel.GetLevelID();
250 std::string fileName = "coordinates_before_rebalance_level_" + buf.str() + ".m";
251 Xpetra::IO<typename Teuchos::ScalarTraits<Scalar>::coordinateType,LO,GO,NO>::Write(fileName,*fineCoords);
252 }
253 if (writeStart <= coarseLevel.GetLevelID() && coarseLevel.GetLevelID() <= writeEnd) {
254 std::ostringstream buf;
255 buf << coarseLevel.GetLevelID();
256 std::string fileName = "coordinates_before_rebalance_level_" + buf.str() + ".m";
257 Xpetra::IO<typename Teuchos::ScalarTraits<Scalar>::coordinateType,LO,GO,NO>::Write(fileName,*coarseCoords);
258 }
259 }
260
261} // namespace MueLu
262
263#endif // MUELU_COORDINATESTRANSFER_FACTORY_DEF_HPP
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
void DeclareInput(Level &finelevel, Level &coarseLevel) const
Specifies the data that this class needs, and the factories that generate that data.
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
Namespace for MueLu classes and methods.
@ Runtime0
One-liner description of what is happening.