MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_GeometricInterpolationPFactory_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_GEOMETRICINTERPOLATIONPFACTORY_DEF_HPP
47 #define MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEF_HPP
48 
49 #include "Xpetra_CrsGraph.hpp"
50 #include "Xpetra_CrsMatrixUtils.hpp"
51 
52 #include "MueLu_MasterList.hpp"
53 #include "MueLu_Monitor.hpp"
54 #include "MueLu_Aggregates.hpp"
55 #include "Xpetra_TpetraCrsMatrix.hpp"
56 
57 // Including this one last ensure that the short names of the above headers are defined properly
58 #include "MueLu_GeometricInterpolationPFactory_decl.hpp"
59 
60 namespace MueLu {
61 
62 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
63 RCP<const ParameterList> GeometricInterpolationPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetValidParameterList() const {
64  RCP<ParameterList> validParamList = rcp(new ParameterList());
65 
66 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
67  SET_VALID_ENTRY("interp: build coarse coordinates");
68 #undef SET_VALID_ENTRY
69 
70  // general variables needed in GeometricInterpolationPFactory
71  validParamList->set<RCP<const FactoryBase> >("A", Teuchos::null,
72  "Generating factory of the matrix A");
73  validParamList->set<RCP<const FactoryBase> >("Aggregates", Teuchos::null,
74  "Aggregates generated by StructuredAggregationFactory used to construct a piece-constant prolongator.");
75  validParamList->set<RCP<const FactoryBase> >("prolongatorGraph", Teuchos::null,
76  "Graph generated by StructuredAggregationFactory used to construct a piece-linear prolongator.");
77  validParamList->set<RCP<const FactoryBase> >("Coordinates", Teuchos::null,
78  "Fine level coordinates used to construct piece-wise linear prolongator and coarse level coordinates.");
79  validParamList->set<RCP<const FactoryBase> >("coarseCoordinatesFineMap", Teuchos::null,
80  "map of the coarse coordinates' GIDs as indexed on the fine mesh.");
81  validParamList->set<RCP<const FactoryBase> >("coarseCoordinatesMap", Teuchos::null,
82  "map of the coarse coordinates' GIDs as indexed on the coarse mesh.");
83  validParamList->set<RCP<const FactoryBase> >("Nullspace", Teuchos::null,
84  "Fine level nullspace used to construct the coarse level nullspace.");
85  validParamList->set<RCP<const FactoryBase> >("numDimensions", Teuchos::null,
86  "Number of spacial dimensions in the problem.");
87  validParamList->set<RCP<const FactoryBase> >("lCoarseNodesPerDim", Teuchos::null,
88  "Number of nodes per spatial dimension on the coarse grid.");
89  validParamList->set<RCP<const FactoryBase> >("structuredInterpolationOrder", Teuchos::null,
90  "Interpolation order for constructing the prolongator.");
91  validParamList->set<bool>("keep coarse coords", false, "Flag to keep coordinates for special coarse grid solve");
92  validParamList->set<bool>("interp: remove small entries", true, "Remove small interpolation coeficient from prolongator to reduce fill-in on coarse level");
93 
94  return validParamList;
95 }
96 
97 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
98 void GeometricInterpolationPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
99  DeclareInput(Level& fineLevel, Level& /* coarseLevel */) const {
100  const ParameterList& pL = GetParameterList();
101 
102  Input(fineLevel, "A");
103  Input(fineLevel, "Nullspace");
104  Input(fineLevel, "numDimensions");
105  Input(fineLevel, "prolongatorGraph");
106  Input(fineLevel, "lCoarseNodesPerDim");
107  Input(fineLevel, "structuredInterpolationOrder");
108 
109  if (pL.get<bool>("interp: build coarse coordinates") ||
110  Get<int>(fineLevel, "structuredInterpolationOrder") == 1) {
111  Input(fineLevel, "Coordinates");
112  Input(fineLevel, "coarseCoordinatesFineMap");
113  Input(fineLevel, "coarseCoordinatesMap");
114  }
115 }
116 
117 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
118 void GeometricInterpolationPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
119  Build(Level& fineLevel, Level& coarseLevel) const {
120  return BuildP(fineLevel, coarseLevel);
121 }
122 
123 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
124 void GeometricInterpolationPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
125  BuildP(Level& fineLevel, Level& coarseLevel) const {
126  FactoryMonitor m(*this, "BuildP", coarseLevel);
127 
128  // Set debug outputs based on environment variable
129  RCP<Teuchos::FancyOStream> out;
130  if (const char* dbg = std::getenv("MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEBUG")) {
131  out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
132  out->setShowAllFrontMatter(false).setShowProcRank(true);
133  } else {
134  out = Teuchos::getFancyOStream(rcp(new Teuchos::oblackholestream()));
135  }
136 
137  *out << "Starting GeometricInterpolationPFactory::BuildP." << std::endl;
138 
139  // Get inputs from the parameter list
140  const ParameterList& pL = GetParameterList();
141  const bool removeSmallEntries = pL.get<bool>("interp: remove small entries");
142  const bool buildCoarseCoordinates = pL.get<bool>("interp: build coarse coordinates");
143  const int interpolationOrder = Get<int>(fineLevel, "structuredInterpolationOrder");
144  const int numDimensions = Get<int>(fineLevel, "numDimensions");
145 
146  // Declared main input/outputs to be retrieved and placed on the fine resp. coarse level
147  RCP<Matrix> A = Get<RCP<Matrix> >(fineLevel, "A");
148  RCP<const CrsGraph> prolongatorGraph = Get<RCP<CrsGraph> >(fineLevel, "prolongatorGraph");
149  RCP<realvaluedmultivector_type> fineCoordinates, coarseCoordinates;
150  RCP<Matrix> P;
151 
152  // Check if we need to build coarse coordinates as they are used if we construct
153  // a linear interpolation prolongator
154  if (buildCoarseCoordinates || (interpolationOrder == 1)) {
155  SubFactoryMonitor sfm(*this, "BuildCoordinates", coarseLevel);
156  RCP<const Map> coarseCoordsFineMap = Get<RCP<const Map> >(fineLevel, "coarseCoordinatesFineMap");
157  RCP<const Map> coarseCoordsMap = Get<RCP<const Map> >(fineLevel, "coarseCoordinatesMap");
158 
159  fineCoordinates = Get<RCP<realvaluedmultivector_type> >(fineLevel, "Coordinates");
160  coarseCoordinates = Xpetra::MultiVectorFactory<real_type, LO, GO, Node>::Build(coarseCoordsFineMap,
161  fineCoordinates->getNumVectors());
162  RCP<const Import> coordsImporter = ImportFactory::Build(fineCoordinates->getMap(),
163  coarseCoordsFineMap);
164  coarseCoordinates->doImport(*fineCoordinates, *coordsImporter, Xpetra::INSERT);
165  coarseCoordinates->replaceMap(coarseCoordsMap);
166 
167  Set(coarseLevel, "Coordinates", coarseCoordinates);
168 
169  if (pL.get<bool>("keep coarse coords")) {
170  coarseLevel.Set<RCP<realvaluedmultivector_type> >("Coordinates2", coarseCoordinates, NoFactory::get());
171  }
172  }
173 
174  *out << "Fine and coarse coordinates have been loaded from the fine level and set on the coarse level." << std::endl;
175 
176  if (interpolationOrder == 0) {
177  SubFactoryMonitor sfm(*this, "BuildConstantP", coarseLevel);
178  // Compute the prolongator using piece-wise constant interpolation
179  BuildConstantP(P, prolongatorGraph, A);
180  } else if (interpolationOrder == 1) {
181  // Compute the prolongator using piece-wise linear interpolation
182  // First get all the required coordinates to compute the local part of P
183  RCP<const Map> prolongatorColMap = prolongatorGraph->getColMap();
184 
185  const size_t dofsPerNode = static_cast<size_t>(A->GetFixedBlockSize());
186  const size_t numColIndices = prolongatorColMap->getLocalNumElements();
187  TEUCHOS_TEST_FOR_EXCEPTION((numColIndices % dofsPerNode) != 0,
188  Exceptions::RuntimeError,
189  "Something went wrong, the number of columns in the prolongator is not a multiple of dofsPerNode!");
190  const size_t numGhostCoords = numColIndices / dofsPerNode;
191  const GO indexBase = prolongatorColMap->getIndexBase();
192  const GO coordIndexBase = fineCoordinates->getMap()->getIndexBase();
193 
194  ArrayView<const GO> prolongatorColIndices = prolongatorColMap->getLocalElementList();
195  Array<GO> ghostCoordIndices(numGhostCoords);
196  for (size_t ghostCoordIdx = 0; ghostCoordIdx < numGhostCoords; ++ghostCoordIdx) {
197  ghostCoordIndices[ghostCoordIdx] = (prolongatorColIndices[ghostCoordIdx * dofsPerNode] - indexBase) / dofsPerNode + coordIndexBase;
198  }
199  RCP<Map> ghostCoordMap = MapFactory::Build(fineCoordinates->getMap()->lib(),
200  prolongatorColMap->getGlobalNumElements() / dofsPerNode,
201  ghostCoordIndices(),
202  coordIndexBase,
203  fineCoordinates->getMap()->getComm());
204 
205  RCP<realvaluedmultivector_type> ghostCoordinates = Xpetra::MultiVectorFactory<real_type, LO, GO, NO>::Build(ghostCoordMap,
206  fineCoordinates->getNumVectors());
207  RCP<const Import> ghostImporter = ImportFactory::Build(coarseCoordinates->getMap(),
208  ghostCoordMap);
209  ghostCoordinates->doImport(*coarseCoordinates, *ghostImporter, Xpetra::INSERT);
210 
211  BuildLinearP(coarseLevel, A, prolongatorGraph, fineCoordinates,
212  ghostCoordinates, numDimensions, removeSmallEntries, P);
213  }
214 
215  *out << "The prolongator matrix has been built." << std::endl;
216 
217  {
218  SubFactoryMonitor sfm(*this, "BuildNullspace", coarseLevel);
219  // Build the coarse nullspace
220  RCP<MultiVector> fineNullspace = Get<RCP<MultiVector> >(fineLevel, "Nullspace");
221  RCP<MultiVector> coarseNullspace = MultiVectorFactory::Build(P->getDomainMap(),
222  fineNullspace->getNumVectors());
223 
224  using helpers = Xpetra::Helpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
225  if (helpers::isTpetraBlockCrs(A)) {
226  // FIXME: BlockCrs doesn't currently support transpose apply, so we have to do this the hard way
227  RCP<Matrix> Ptrans = Utilities::Transpose(*P);
228  Ptrans->apply(*fineNullspace, *coarseNullspace, Teuchos::NO_TRANS, Teuchos::ScalarTraits<SC>::one(),
229  Teuchos::ScalarTraits<SC>::zero());
230  } else {
231  P->apply(*fineNullspace, *coarseNullspace, Teuchos::TRANS, Teuchos::ScalarTraits<SC>::one(),
232  Teuchos::ScalarTraits<SC>::zero());
233  }
234 
235  Set(coarseLevel, "Nullspace", coarseNullspace);
236  }
237 
238  *out << "The coarse nullspace is constructed and set on the coarse level." << std::endl;
239 
240  Set(coarseLevel, "P", P);
241 
242  *out << "GeometricInterpolationPFactory::BuildP has completed." << std::endl;
243 
244 } // BuildP
245 
246 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
247 void GeometricInterpolationPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
248  BuildConstantP(RCP<Matrix>& P, RCP<const CrsGraph>& prolongatorGraph, RCP<Matrix>& A) const {
249  // Set debug outputs based on environment variable
250  RCP<Teuchos::FancyOStream> out;
251  if (const char* dbg = std::getenv("MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEBUG")) {
252  out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
253  out->setShowAllFrontMatter(false).setShowProcRank(true);
254  } else {
255  out = Teuchos::getFancyOStream(rcp(new Teuchos::oblackholestream()));
256  }
257 
258  *out << "BuildConstantP" << std::endl;
259 
260  std::vector<size_t> strideInfo(1);
261  strideInfo[0] = A->GetFixedBlockSize();
262  RCP<const StridedMap> stridedDomainMap =
263  StridedMapFactory::Build(prolongatorGraph->getDomainMap(), strideInfo);
264 
265  *out << "Call prolongator constructor" << std::endl;
266 
267  using helpers = Xpetra::Helpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
268  if (helpers::isTpetraBlockCrs(A)) {
269  SC one = Teuchos::ScalarTraits<SC>::one();
270  SC zero = Teuchos::ScalarTraits<SC>::zero();
271  LO NSDim = A->GetStorageBlockSize();
272 
273  // Build the exploded Map
274  RCP<const Map> BlockMap = prolongatorGraph->getDomainMap();
275  Teuchos::ArrayView<const GO> block_dofs = BlockMap->getLocalElementList();
276  Teuchos::Array<GO> point_dofs(block_dofs.size() * NSDim);
277  for (LO i = 0, ct = 0; i < block_dofs.size(); i++) {
278  for (LO j = 0; j < NSDim; j++) {
279  point_dofs[ct] = block_dofs[i] * NSDim + j;
280  ct++;
281  }
282  }
283 
284  RCP<const Map> PointMap = MapFactory::Build(BlockMap->lib(),
285  BlockMap->getGlobalNumElements() * NSDim,
286  point_dofs(),
287  BlockMap->getIndexBase(),
288  BlockMap->getComm());
289  strideInfo[0] = A->GetFixedBlockSize();
290  RCP<const StridedMap> stridedPointMap = StridedMapFactory::Build(PointMap, strideInfo);
291 
292  RCP<Xpetra::CrsMatrix<SC, LO, GO, NO> > P_xpetra = Xpetra::CrsMatrixFactory<SC, LO, GO, NO>::BuildBlock(prolongatorGraph, PointMap, A->getRangeMap(), NSDim);
293  RCP<Xpetra::TpetraBlockCrsMatrix<SC, LO, GO, NO> > P_tpetra = rcp_dynamic_cast<Xpetra::TpetraBlockCrsMatrix<SC, LO, GO, NO> >(P_xpetra);
294  if (P_tpetra.is_null()) throw std::runtime_error("BuildConstantP Matrix factory did not return a Tpetra::BlockCrsMatrix");
295  RCP<CrsMatrixWrap> P_wrap = rcp(new CrsMatrixWrap(P_xpetra));
296 
297  // NOTE: Assumes block-diagonal prolongation
298  Teuchos::Array<LO> temp(1);
299  Teuchos::ArrayView<const LO> indices;
300  Teuchos::Array<Scalar> block(NSDim * NSDim, zero);
301  for (LO i = 0; i < NSDim; i++)
302  block[i * NSDim + i] = one;
303  for (LO i = 0; i < (int)prolongatorGraph->getLocalNumRows(); i++) {
304  prolongatorGraph->getLocalRowView(i, indices);
305  for (LO j = 0; j < (LO)indices.size(); j++) {
306  temp[0] = indices[j];
307  P_tpetra->replaceLocalValues(i, temp(), block());
308  }
309  }
310 
311  P = P_wrap;
312  if (A->IsView("stridedMaps") == true) {
313  P->CreateView("stridedMaps", A->getRowMap("stridedMaps"), stridedPointMap);
314  } else {
315  P->CreateView("stridedMaps", P->getRangeMap(), PointMap);
316  }
317  } else {
318  // Create the prolongator matrix and its associated objects
319  RCP<ParameterList> dummyList = rcp(new ParameterList());
320  P = rcp(new CrsMatrixWrap(prolongatorGraph, dummyList));
321  RCP<CrsMatrix> PCrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
322  PCrs->setAllToScalar(1.0);
323  PCrs->fillComplete();
324 
325  // set StridingInformation of P
326  if (A->IsView("stridedMaps") == true)
327  P->CreateView("stridedMaps", A->getRowMap("stridedMaps"), stridedDomainMap);
328  else
329  P->CreateView("stridedMaps", P->getRangeMap(), stridedDomainMap);
330  }
331 
332 } // BuildConstantP
333 
334 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
335 void GeometricInterpolationPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
336  BuildLinearP(Level& coarseLevel,
337  RCP<Matrix>& A, RCP<const CrsGraph>& prolongatorGraph,
338  RCP<realvaluedmultivector_type>& fineCoordinates,
339  RCP<realvaluedmultivector_type>& ghostCoordinates,
340  const int numDimensions, const bool removeSmallEntries,
341  RCP<Matrix>& P) const {
342  // Set debug outputs based on environment variable
343  RCP<Teuchos::FancyOStream> out;
344  if (const char* dbg = std::getenv("MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEBUG")) {
345  out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
346  out->setShowAllFrontMatter(false).setShowProcRank(true);
347  } else {
348  out = Teuchos::getFancyOStream(rcp(new Teuchos::oblackholestream()));
349  }
350 
351  // Compute 2^numDimensions using bit logic to avoid round-off errors
352  const int numInterpolationPoints = 1 << numDimensions;
353  const int dofsPerNode = A->GetFixedBlockSize() / A->GetStorageBlockSize();
354  ;
355 
356  RCP<ParameterList> dummyList = rcp(new ParameterList());
357  P = rcp(new CrsMatrixWrap(prolongatorGraph, dummyList));
358  RCP<CrsMatrix> PCrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
359  PCrs->resumeFill(); // The Epetra matrix is considered filled at this point.
360 
361  {
362  *out << "Entering BuildLinearP" << std::endl;
363  SubFactoryMonitor sfm(*this, "BuildLinearP", coarseLevel);
364 
365  // Extract coordinates for interpolation stencil calculations
366  const LO numFineNodes = fineCoordinates->getLocalLength();
367  const LO numGhostNodes = ghostCoordinates->getLocalLength();
368  Array<ArrayRCP<const real_type> > fineCoords(3);
369  Array<ArrayRCP<const real_type> > ghostCoords(3);
370  const real_type realZero = Teuchos::as<real_type>(0.0);
371  ArrayRCP<real_type> fineZero(numFineNodes, realZero);
372  ArrayRCP<real_type> ghostZero(numGhostNodes, realZero);
373  for (int dim = 0; dim < 3; ++dim) {
374  if (dim < numDimensions) {
375  fineCoords[dim] = fineCoordinates->getData(dim);
376  ghostCoords[dim] = ghostCoordinates->getData(dim);
377  } else {
378  fineCoords[dim] = fineZero;
379  ghostCoords[dim] = ghostZero;
380  }
381  }
382 
383  *out << "Coordinates extracted from the multivectors!" << std::endl;
384 
385  { // Construct the linear interpolation prolongator
386  LO interpolationNodeIdx = 0, rowIdx = 0;
387  ArrayView<const LO> colIndices;
388  Array<SC> values;
389  Array<Array<real_type> > coords(numInterpolationPoints + 1);
390  Array<real_type> stencil(numInterpolationPoints);
391  for (LO nodeIdx = 0; nodeIdx < numFineNodes; ++nodeIdx) {
392  if (PCrs->getNumEntriesInLocalRow(nodeIdx * dofsPerNode) == 1) {
393  values.resize(1);
394  values[0] = 1.0;
395  for (LO dof = 0; dof < dofsPerNode; ++dof) {
396  rowIdx = nodeIdx * dofsPerNode + dof;
397  prolongatorGraph->getLocalRowView(rowIdx, colIndices);
398  PCrs->replaceLocalValues(rowIdx, colIndices, values());
399  }
400  } else {
401  // Extract the coordinates associated with the current node
402  // and the neighboring coarse nodes
403  coords[0].resize(3);
404  for (int dim = 0; dim < 3; ++dim) {
405  coords[0][dim] = fineCoords[dim][nodeIdx];
406  }
407  prolongatorGraph->getLocalRowView(nodeIdx * dofsPerNode, colIndices);
408  for (int interpolationIdx = 0; interpolationIdx < numInterpolationPoints; ++interpolationIdx) {
409  coords[interpolationIdx + 1].resize(3);
410  interpolationNodeIdx = colIndices[interpolationIdx] / dofsPerNode;
411  for (int dim = 0; dim < 3; ++dim) {
412  coords[interpolationIdx + 1][dim] = ghostCoords[dim][interpolationNodeIdx];
413  }
414  }
415  RCP<Teuchos::TimeMonitor> tm = rcp(new Teuchos::TimeMonitor(*Teuchos::TimeMonitor::getNewTimer("Compute Linear Interpolation")));
416  ComputeLinearInterpolationStencil(numDimensions, numInterpolationPoints, coords, stencil);
417  tm = Teuchos::null;
418  values.resize(numInterpolationPoints);
419  for (LO valueIdx = 0; valueIdx < numInterpolationPoints; ++valueIdx) {
420  values[valueIdx] = Teuchos::as<SC>(stencil[valueIdx]);
421  }
422 
423  // Set values in all the rows corresponding to nodeIdx
424  for (LO dof = 0; dof < dofsPerNode; ++dof) {
425  rowIdx = nodeIdx * dofsPerNode + dof;
426  prolongatorGraph->getLocalRowView(rowIdx, colIndices);
427  PCrs->replaceLocalValues(rowIdx, colIndices, values());
428  }
429  } // Check for Coarse vs. Fine point
430  } // Loop over fine nodes
431  } // Construct the linear interpolation prolongator
432 
433  *out << "The calculation of the interpolation stencils has completed." << std::endl;
434 
435  PCrs->fillComplete();
436  }
437 
438  *out << "All values in P have been set and fillComplete has been performed." << std::endl;
439 
440  // Note lbv Jan 29 2019: this should be handle at aggregation level
441  // if the user really does not want potential d2 neighbors on coarse grid
442  // that way we would avoid a new graph construction...
443 
444  // Check if we want to remove small entries from P
445  // to reduce stencil growth on next level.
446  if (removeSmallEntries) {
447  *out << "Entering remove small entries" << std::endl;
448  SubFactoryMonitor sfm(*this, "remove small entries", coarseLevel);
449 
450  ArrayRCP<const size_t> rowptrOrig;
451  ArrayRCP<const LO> colindOrig;
452  ArrayRCP<const Scalar> valuesOrig;
453  PCrs->getAllValues(rowptrOrig, colindOrig, valuesOrig);
454 
455  const size_t numRows = static_cast<size_t>(rowptrOrig.size() - 1);
456  ArrayRCP<size_t> rowPtr(numRows + 1);
457  ArrayRCP<size_t> nnzOnRows(numRows);
458  rowPtr[0] = 0;
459  size_t countRemovedEntries = 0;
460  for (size_t rowIdx = 0; rowIdx < numRows; ++rowIdx) {
461  for (size_t entryIdx = rowptrOrig[rowIdx]; entryIdx < rowptrOrig[rowIdx + 1]; ++entryIdx) {
462  if (Teuchos::ScalarTraits<Scalar>::magnitude(valuesOrig[entryIdx]) < 1e-6) {
463  ++countRemovedEntries;
464  }
465  }
466  rowPtr[rowIdx + 1] = rowptrOrig[rowIdx + 1] - countRemovedEntries;
467  nnzOnRows[rowIdx] = rowPtr[rowIdx + 1] - rowPtr[rowIdx];
468  }
469  GetOStream(Statistics1) << "interp: number of small entries removed= " << countRemovedEntries << " / " << rowptrOrig[numRows] << std::endl;
470 
471  size_t countKeptEntries = 0;
472  ArrayRCP<LO> colInd(rowPtr[numRows]);
473  ArrayRCP<SC> values(rowPtr[numRows]);
474  for (size_t entryIdx = 0; entryIdx < rowptrOrig[numRows]; ++entryIdx) {
475  if (Teuchos::ScalarTraits<Scalar>::magnitude(valuesOrig[entryIdx]) > 1e-6) {
476  colInd[countKeptEntries] = colindOrig[entryIdx];
477  values[countKeptEntries] = valuesOrig[entryIdx];
478  ++countKeptEntries;
479  }
480  }
481 
482  P = rcp(new CrsMatrixWrap(prolongatorGraph->getRowMap(),
483  prolongatorGraph->getColMap(),
484  nnzOnRows));
485  RCP<CrsMatrix> PCrsSqueezed = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
486  PCrsSqueezed->resumeFill(); // The Epetra matrix is considered filled at this point.
487  PCrsSqueezed->setAllValues(rowPtr, colInd, values);
488  PCrsSqueezed->expertStaticFillComplete(prolongatorGraph->getDomainMap(),
489  prolongatorGraph->getRangeMap());
490  }
491 
492  std::vector<size_t> strideInfo(1);
493  strideInfo[0] = dofsPerNode;
494  RCP<const StridedMap> stridedDomainMap =
495  StridedMapFactory::Build(prolongatorGraph->getDomainMap(), strideInfo);
496 
497  *out << "The strided maps of P have been computed" << std::endl;
498 
499  // set StridingInformation of P
500  if (A->IsView("stridedMaps") == true) {
501  P->CreateView("stridedMaps", A->getRowMap("stridedMaps"), stridedDomainMap);
502  } else {
503  P->CreateView("stridedMaps", P->getRangeMap(), stridedDomainMap);
504  }
505 
506 } // BuildLinearP
507 
508 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
509 void GeometricInterpolationPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
510  ComputeLinearInterpolationStencil(const int numDimensions, const int numInterpolationPoints,
511  const Array<Array<real_type> > coord,
512  Array<real_type>& stencil) const {
513  // 7 8 Find xi, eta and zeta such that
514  // x---------x
515  // /| /| Rx = x_p - sum N_i(xi,eta,zeta)x_i = 0
516  // 5/ | 6/ | Ry = y_p - sum N_i(xi,eta,zeta)y_i = 0
517  // x---------x | Rz = z_p - sum N_i(xi,eta,zeta)z_i = 0
518  // | | *P | |
519  // | x------|--x We can do this with a Newton solver:
520  // | /3 | /4 We will start with initial guess (xi,eta,zeta) = (0,0,0)
521  // |/ |/ We compute the Jacobian and iterate until convergence...
522  // z y x---------x
523  // | / 1 2 Once we have (xi,eta,zeta), we can evaluate all N_i which
524  // |/ give us the weights for the interpolation stencil!
525  // o---x
526  //
527 
528  Teuchos::SerialDenseMatrix<LO, real_type> Jacobian(numDimensions, numDimensions);
529  Teuchos::SerialDenseVector<LO, real_type> residual(numDimensions);
530  Teuchos::SerialDenseVector<LO, real_type> solutionDirection(numDimensions);
531  Teuchos::SerialDenseVector<LO, real_type> paramCoords(numDimensions);
532  Teuchos::SerialDenseSolver<LO, real_type> problem;
533  int iter = 0, max_iter = 5;
534  real_type functions[4][8], norm_ref = 1.0, norm2 = 1.0, tol = 1.0e-5;
535  paramCoords.size(numDimensions);
536 
537  while ((iter < max_iter) && (norm2 > tol * norm_ref)) {
538  ++iter;
539  norm2 = 0.0;
540  solutionDirection.size(numDimensions);
541  residual.size(numDimensions);
542  Jacobian = 0.0;
543 
544  // Compute Jacobian and Residual
545  GetInterpolationFunctions(numDimensions, paramCoords, functions);
546  for (LO i = 0; i < numDimensions; ++i) {
547  residual(i) = coord[0][i]; // Add coordinates from point of interest
548  for (LO k = 0; k < numInterpolationPoints; ++k) {
549  residual(i) -= functions[0][k] * coord[k + 1][i]; // Remove contribution from all coarse points
550  }
551  if (iter == 1) {
552  norm_ref += residual(i) * residual(i);
553  if (i == numDimensions - 1) {
554  norm_ref = std::sqrt(norm_ref);
555  }
556  }
557 
558  for (LO j = 0; j < numDimensions; ++j) {
559  for (LO k = 0; k < numInterpolationPoints; ++k) {
560  Jacobian(i, j) += functions[j + 1][k] * coord[k + 1][i];
561  }
562  }
563  }
564 
565  // Set Jacobian, Vectors and solve problem
566  problem.setMatrix(Teuchos::rcp(&Jacobian, false));
567  problem.setVectors(Teuchos::rcp(&solutionDirection, false), Teuchos::rcp(&residual, false));
568  if (problem.shouldEquilibrate()) {
569  problem.factorWithEquilibration(true);
570  }
571  problem.solve();
572 
573  for (LO i = 0; i < numDimensions; ++i) {
574  paramCoords(i) = paramCoords(i) + solutionDirection(i);
575  }
576 
577  // Recompute Residual norm
578  GetInterpolationFunctions(numDimensions, paramCoords, functions);
579  for (LO i = 0; i < numDimensions; ++i) {
580  real_type tmp = coord[0][i];
581  for (LO k = 0; k < numInterpolationPoints; ++k) {
582  tmp -= functions[0][k] * coord[k + 1][i];
583  }
584  norm2 += tmp * tmp;
585  tmp = 0.0;
586  }
587  norm2 = std::sqrt(norm2);
588  }
589 
590  // Load the interpolation values onto the stencil.
591  for (LO i = 0; i < numInterpolationPoints; ++i) {
592  stencil[i] = functions[0][i];
593  }
594 
595 } // End ComputeLinearInterpolationStencil
596 
597 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
598 void GeometricInterpolationPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
599  GetInterpolationFunctions(const LO numDimensions,
600  const Teuchos::SerialDenseVector<LO, real_type> parametricCoordinates,
601  real_type functions[4][8]) const {
602  real_type xi = 0.0, eta = 0.0, zeta = 0.0, denominator = 0.0;
603  if (numDimensions == 1) {
604  xi = parametricCoordinates[0];
605  denominator = 2.0;
606  } else if (numDimensions == 2) {
607  xi = parametricCoordinates[0];
608  eta = parametricCoordinates[1];
609  denominator = 4.0;
610  } else if (numDimensions == 3) {
611  xi = parametricCoordinates[0];
612  eta = parametricCoordinates[1];
613  zeta = parametricCoordinates[2];
614  denominator = 8.0;
615  }
616 
617  functions[0][0] = (1.0 - xi) * (1.0 - eta) * (1.0 - zeta) / denominator;
618  functions[0][1] = (1.0 + xi) * (1.0 - eta) * (1.0 - zeta) / denominator;
619  functions[0][2] = (1.0 - xi) * (1.0 + eta) * (1.0 - zeta) / denominator;
620  functions[0][3] = (1.0 + xi) * (1.0 + eta) * (1.0 - zeta) / denominator;
621  functions[0][4] = (1.0 - xi) * (1.0 - eta) * (1.0 + zeta) / denominator;
622  functions[0][5] = (1.0 + xi) * (1.0 - eta) * (1.0 + zeta) / denominator;
623  functions[0][6] = (1.0 - xi) * (1.0 + eta) * (1.0 + zeta) / denominator;
624  functions[0][7] = (1.0 + xi) * (1.0 + eta) * (1.0 + zeta) / denominator;
625 
626  functions[1][0] = -(1.0 - eta) * (1.0 - zeta) / denominator;
627  functions[1][1] = (1.0 - eta) * (1.0 - zeta) / denominator;
628  functions[1][2] = -(1.0 + eta) * (1.0 - zeta) / denominator;
629  functions[1][3] = (1.0 + eta) * (1.0 - zeta) / denominator;
630  functions[1][4] = -(1.0 - eta) * (1.0 + zeta) / denominator;
631  functions[1][5] = (1.0 - eta) * (1.0 + zeta) / denominator;
632  functions[1][6] = -(1.0 + eta) * (1.0 + zeta) / denominator;
633  functions[1][7] = (1.0 + eta) * (1.0 + zeta) / denominator;
634 
635  functions[2][0] = -(1.0 - xi) * (1.0 - zeta) / denominator;
636  functions[2][1] = -(1.0 + xi) * (1.0 - zeta) / denominator;
637  functions[2][2] = (1.0 - xi) * (1.0 - zeta) / denominator;
638  functions[2][3] = (1.0 + xi) * (1.0 - zeta) / denominator;
639  functions[2][4] = -(1.0 - xi) * (1.0 + zeta) / denominator;
640  functions[2][5] = -(1.0 + xi) * (1.0 + zeta) / denominator;
641  functions[2][6] = (1.0 - xi) * (1.0 + zeta) / denominator;
642  functions[2][7] = (1.0 + xi) * (1.0 + zeta) / denominator;
643 
644  functions[3][0] = -(1.0 - xi) * (1.0 - eta) / denominator;
645  functions[3][1] = -(1.0 + xi) * (1.0 - eta) / denominator;
646  functions[3][2] = -(1.0 - xi) * (1.0 + eta) / denominator;
647  functions[3][3] = -(1.0 + xi) * (1.0 + eta) / denominator;
648  functions[3][4] = (1.0 - xi) * (1.0 - eta) / denominator;
649  functions[3][5] = (1.0 + xi) * (1.0 - eta) / denominator;
650  functions[3][6] = (1.0 - xi) * (1.0 + eta) / denominator;
651  functions[3][7] = (1.0 + xi) * (1.0 + eta) / denominator;
652 
653 } // End GetInterpolationFunctions
654 
655 } // namespace MueLu
656 
657 #endif // MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEF_HPP
Xpetra::MultiVectorFactory::Build
static Teuchos::RCP< MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Build(const Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node >> &map, size_t NumVectors, bool zeroOut=true)
Teuchos::basic_FancyOStream::setShowProcRank
basic_FancyOStream & setShowProcRank(const bool showProcRank)
GO
GlobalOrdinal GO
Teuchos::ParameterList::get
T & get(const std::string &name, T def_value)
Teuchos::ParameterList::set
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
MueLu::FactoryMonitor
Timer to be used in factories. Similar to Monitor but with additional timers.
Definition: MueLu_Monitor.hpp:199
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Teuchos::basic_oblackholestream
MueLu::Statistics1
Print more statistics.
Definition: MueLu_VerbosityLevel.hpp:70
Teuchos::ArrayView::size
size_type size() const
MueLu::GeometricInterpolationPFactory::BuildLinearP
void BuildLinearP(Level &coarseLevel, RCP< Matrix > &A, RCP< const CrsGraph > &prolongatorGraph, RCP< realvaluedmultivector_type > &fineCoordinates, RCP< realvaluedmultivector_type > &ghostCoordinates, const int numDimensions, const bool keepD2, RCP< Matrix > &P) const
Definition: MueLu_GeometricInterpolationPFactory_def.hpp:336
LO
LocalOrdinal LO
MueLu_GeometricInterpolationPFactory_decl.hpp
Teuchos::ArrayRCP::size
size_type size() const
MueLu::NoFactory::get
static const NoFactory * get()
Definition: MueLu_NoFactory.cpp:59
MueLu::GeometricInterpolationPFactory::GetInterpolationFunctions
void GetInterpolationFunctions(const LO numDimensions, const Teuchos::SerialDenseVector< LO, real_type > parametricCoordinates, real_type functions[4][8]) const
Definition: MueLu_GeometricInterpolationPFactory_def.hpp:599
Xpetra_CrsGraph.hpp
MueLu_MasterList.hpp
Teuchos::TimeMonitor::getNewTimer
static RCP< Time > getNewTimer(const std::string &name)
rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
MueLu::GeometricInterpolationPFactory::Build
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
Definition: MueLu_GeometricInterpolationPFactory_def.hpp:119
MueLu::Level
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
MueLu::GeometricInterpolationPFactory::ComputeLinearInterpolationStencil
void ComputeLinearInterpolationStencil(const int numDimensions, const int numInterpolationPoints, const Array< Array< real_type > > coord, Array< real_type > &stencil) const
Definition: MueLu_GeometricInterpolationPFactory_def.hpp:510
MueLu::SubFactoryMonitor
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
Definition: MueLu_Monitor.hpp:250
Teuchos::SerialDenseSolver::factorWithEquilibration
void factorWithEquilibration(bool flag)
SET_VALID_ENTRY
#define SET_VALID_ENTRY(name)
MueLu::Utilities::Transpose
static RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Transpose(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, bool optimizeTranspose=false, const std::string &label=std::string(), const Teuchos::RCP< Teuchos::ParameterList > &params=Teuchos::null)
Definition: MueLu_Utilities_def.hpp:532
Teuchos::basic_FancyOStream::setShowAllFrontMatter
basic_FancyOStream & setShowAllFrontMatter(const bool showAllFrontMatter)
Teuchos::Array::resize
void resize(size_type new_size, const value_type &x=value_type())
MueLu::Level::Set
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
Definition: MueLu_Level.hpp:161
MueLu::GeometricInterpolationPFactory::DeclareInput
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
Definition: MueLu_GeometricInterpolationPFactory_def.hpp:99
iter
TransListIter iter
Teuchos::SerialDenseVector::size
int size(OrdinalType length_in)
Teuchos::SerialDenseSolver::setVectors
int setVectors(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &X, const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &B)
MueLu::GeometricInterpolationPFactory::GetValidParameterList
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Definition: MueLu_GeometricInterpolationPFactory_def.hpp:63
SC
Scalar SC
MueLu_Monitor.hpp
residual
void residual(const Operator< SC, LO, GO, NO > &Aop, const MultiVector< SC, LO, GO, NO > &X_in, const MultiVector< SC, LO, GO, NO > &B_in, MultiVector< SC, LO, GO, NO > &R_in)
MueLu::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition: MueLu_Exceptions.hpp:70
MueLu::GeometricInterpolationPFactory::BuildConstantP
void BuildConstantP(RCP< Matrix > &P, RCP< const CrsGraph > &prolongatorGraph, RCP< Matrix > &A) const
Definition: MueLu_GeometricInterpolationPFactory_def.hpp:248
MueLu::GeometricInterpolationPFactory::real_type
typename Teuchos::ScalarTraits< SC >::coordinateType real_type
Definition: MueLu_GeometricInterpolationPFactory_decl.hpp:69
Xpetra_CrsMatrixUtils.hpp
MueLu::GeometricInterpolationPFactory::BuildP
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
Definition: MueLu_GeometricInterpolationPFactory_def.hpp:125
Teuchos::SerialDenseSolver::setMatrix
int setMatrix(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &A)