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