46 #ifndef MUELU_GEOMETRICINTERPOLATIONPFACTORY_KOKKOS_DEF_HPP 47 #define MUELU_GEOMETRICINTERPOLATIONPFACTORY_KOKKOS_DEF_HPP 49 #include "Xpetra_CrsGraph.hpp" 50 #include "Xpetra_CrsMatrixUtils.hpp" 54 #include "MueLu_IndexManager_kokkos.hpp" 61 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
63 RCP<ParameterList> validParamList = rcp(
new ParameterList());
65 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name)) 67 #undef SET_VALID_ENTRY 70 validParamList->set<RCP<const FactoryBase> >(
"A", Teuchos::null,
71 "Generating factory of the matrix A");
72 validParamList->set<RCP<const FactoryBase> >(
"prolongatorGraph", Teuchos::null,
73 "Graph generated by StructuredAggregationFactory used to construct a piece-linear prolongator.");
74 validParamList->set<RCP<const FactoryBase> >(
"Coordinates", Teuchos::null,
75 "Fine level coordinates used to construct piece-wise linear prolongator and coarse level coordinates.");
76 validParamList->set<RCP<const FactoryBase> >(
"Nullspace", Teuchos::null,
77 "Fine level nullspace used to construct the coarse level nullspace.");
78 validParamList->set<RCP<const FactoryBase> >(
"numDimensions", Teuchos::null,
79 "Number of spacial dimensions in the problem.");
80 validParamList->set<RCP<const FactoryBase> >(
"lCoarseNodesPerDim", Teuchos::null,
81 "Number of nodes per spatial dimension on the coarse grid.");
82 validParamList->set<RCP<const FactoryBase> >(
"indexManager", Teuchos::null,
83 "The index manager associated with the local mesh.");
84 validParamList->set<RCP<const FactoryBase> >(
"structuredInterpolationOrder", Teuchos::null,
85 "Interpolation order for constructing the prolongator.");
87 return validParamList;
90 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
93 const ParameterList& pL = GetParameterList();
95 Input(fineLevel,
"A");
96 Input(fineLevel,
"Nullspace");
97 Input(fineLevel,
"numDimensions");
98 Input(fineLevel,
"prolongatorGraph");
99 Input(fineLevel,
"lCoarseNodesPerDim");
100 Input(fineLevel,
"structuredInterpolationOrder");
102 if( pL.get<
bool>(
"interp: build coarse coordinates") ||
103 Get<int>(fineLevel,
"structuredInterpolationOrder") == 1) {
104 Input(fineLevel,
"Coordinates");
105 Input(fineLevel,
"indexManager");
110 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
113 return BuildP(fineLevel, coarseLevel);
116 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
122 RCP<Teuchos::FancyOStream> out;
123 if(
const char* dbg = std::getenv(
"MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEBUG")) {
124 out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
125 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
127 out = Teuchos::getFancyOStream(rcp(
new Teuchos::oblackholestream()));
130 *out <<
"Starting GeometricInterpolationPFactory_kokkos::BuildP." << std::endl;
133 const ParameterList& pL = GetParameterList();
134 const bool buildCoarseCoordinates = pL.get<
bool>(
"interp: build coarse coordinates");
135 const int interpolationOrder = Get<int>(fineLevel,
"structuredInterpolationOrder");
136 const int numDimensions = Get<int>(fineLevel,
"numDimensions");
139 RCP<Matrix> A = Get<RCP<Matrix> >(fineLevel,
"A");
140 RCP<const CrsGraph> prolongatorGraph = Get<RCP<CrsGraph> >(fineLevel,
"prolongatorGraph");
141 RCP<realvaluedmultivector_type> fineCoordinates, coarseCoordinates;
146 if(buildCoarseCoordinates || (interpolationOrder == 1)) {
150 RCP<IndexManager_kokkos> geoData = Get<RCP<IndexManager_kokkos> >(fineLevel,
"indexManager");
151 fineCoordinates = Get< RCP<realvaluedmultivector_type> >(fineLevel,
"Coordinates");
154 RCP<const Map> coarseCoordsMap = MapFactory::Build(fineCoordinates->getMap()->lib(),
155 Teuchos::OrdinalTraits<GO>::invalid(),
156 geoData->getNumCoarseNodes(),
157 fineCoordinates->getMap()->getIndexBase(),
158 fineCoordinates->getMap()->getComm());
159 coarseCoordinates = Xpetra::MultiVectorFactory<real_type,LO,GO,Node>::
160 Build(coarseCoordsMap, fineCoordinates->getNumVectors());
165 fineCoordinates-> getDeviceLocalView(Xpetra::Access::ReadWrite),
166 coarseCoordinates->getDeviceLocalView(Xpetra::Access::OverwriteAll));
167 Kokkos::parallel_for(
"GeometricInterpolation: build coarse coordinates",
169 myCoarseCoordinatesBuilder);
171 Set(coarseLevel,
"Coordinates", coarseCoordinates);
174 *out <<
"Fine and coarse coordinates have been loaded from the fine level and set on the coarse level." << std::endl;
176 if(interpolationOrder == 0) {
179 BuildConstantP(P, prolongatorGraph, A);
180 }
else if(interpolationOrder == 1) {
183 RCP<realvaluedmultivector_type> ghostCoordinates
184 = Xpetra::MultiVectorFactory<real_type,LO,GO,NO>::Build(prolongatorGraph->getColMap(),
185 fineCoordinates->getNumVectors());
186 RCP<const Import> ghostImporter = ImportFactory::Build(coarseCoordinates->getMap(),
187 prolongatorGraph->getColMap());
188 ghostCoordinates->doImport(*coarseCoordinates, *ghostImporter, Xpetra::INSERT);
191 BuildLinearP(A, prolongatorGraph, fineCoordinates, ghostCoordinates, numDimensions, P);
194 *out <<
"The prolongator matrix has been built." << std::endl;
199 RCP<MultiVector> fineNullspace = Get< RCP<MultiVector> > (fineLevel,
"Nullspace");
200 RCP<MultiVector> coarseNullspace = MultiVectorFactory::Build(P->getDomainMap(),
201 fineNullspace->getNumVectors());
202 P->apply(*fineNullspace, *coarseNullspace, Teuchos::TRANS, Teuchos::ScalarTraits<SC>::one(),
203 Teuchos::ScalarTraits<SC>::zero());
204 Set(coarseLevel,
"Nullspace", coarseNullspace);
207 *out <<
"The coarse nullspace is constructed and set on the coarse level." << std::endl;
209 Array<LO> lNodesPerDir = Get<Array<LO> >(fineLevel,
"lCoarseNodesPerDim");
210 Set(coarseLevel,
"numDimensions", numDimensions);
211 Set(coarseLevel,
"lNodesPerDim", lNodesPerDir);
212 Set(coarseLevel,
"P", P);
214 *out <<
"GeometricInterpolationPFactory_kokkos::BuildP has completed." << std::endl;
218 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
220 BuildConstantP(RCP<Matrix>& P, RCP<const CrsGraph>& prolongatorGraph, RCP<Matrix>& A)
const {
223 RCP<Teuchos::FancyOStream> out;
224 if(
const char* dbg = std::getenv(
"MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEBUG")) {
225 out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
226 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
228 out = Teuchos::getFancyOStream(rcp(
new Teuchos::oblackholestream()));
231 *out <<
"BuildConstantP" << std::endl;
233 std::vector<size_t> strideInfo(1);
234 strideInfo[0] = A->GetFixedBlockSize();
235 RCP<const StridedMap> stridedDomainMap =
236 StridedMapFactory::Build(prolongatorGraph->getDomainMap(), strideInfo);
238 *out <<
"Call prolongator constructor" << std::endl;
241 RCP<ParameterList> dummyList = rcp(
new ParameterList());
242 P = rcp(
new CrsMatrixWrap(prolongatorGraph, dummyList));
243 RCP<CrsMatrix> PCrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
244 PCrs->setAllToScalar(1.0);
245 PCrs->fillComplete();
248 if (A->IsView(
"stridedMaps") ==
true) {
249 P->CreateView(
"stridedMaps", A->getRowMap(
"stridedMaps"), stridedDomainMap);
251 P->CreateView(
"stridedMaps", P->getRangeMap(), stridedDomainMap);
256 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
259 RCP<realvaluedmultivector_type>& fineCoordinates,
260 RCP<realvaluedmultivector_type>& ghostCoordinates,
261 const int numDimensions, RCP<Matrix>& P)
const {
264 RCP<Teuchos::FancyOStream> out;
265 if(
const char* dbg = std::getenv(
"MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEBUG")) {
266 out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
267 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
269 out = Teuchos::getFancyOStream(rcp(
new Teuchos::oblackholestream()));
272 *out <<
"Entering BuildLinearP" << std::endl;
275 const LO numFineNodes = fineCoordinates->getLocalLength();
276 const LO numGhostNodes = ghostCoordinates->getLocalLength();
277 Array<ArrayRCP<const real_type> > fineCoords(3);
278 Array<ArrayRCP<const real_type> > ghostCoords(3);
279 const real_type realZero = Teuchos::as<real_type>(0.0);
280 ArrayRCP<real_type> fineZero(numFineNodes, realZero);
281 ArrayRCP<real_type> ghostZero(numGhostNodes, realZero);
282 for(
int dim = 0; dim < 3; ++dim) {
283 if(dim < numDimensions) {
284 fineCoords[dim] = fineCoordinates->getData(dim);
285 ghostCoords[dim] = ghostCoordinates->getData(dim);
287 fineCoords[dim] = fineZero;
288 ghostCoords[dim] = ghostZero;
292 *out <<
"Coordinates extracted from the multivectors!" << std::endl;
295 const int numInterpolationPoints = 1 << numDimensions;
296 const int dofsPerNode = A->GetFixedBlockSize();
298 std::vector<size_t> strideInfo(1);
299 strideInfo[0] = dofsPerNode;
300 RCP<const StridedMap> stridedDomainMap =
301 StridedMapFactory::Build(prolongatorGraph->getDomainMap(), strideInfo);
303 *out <<
"The maps of P have been computed" << std::endl;
305 RCP<ParameterList> dummyList = rcp(
new ParameterList());
306 P = rcp(
new CrsMatrixWrap(prolongatorGraph, dummyList));
307 RCP<CrsMatrix> PCrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
310 LO interpolationNodeIdx = 0, rowIdx = 0;
311 ArrayView<const LO> colIndices;
313 Array<Array<real_type> > coords(numInterpolationPoints + 1);
314 Array<real_type> stencil(numInterpolationPoints);
315 for(LO nodeIdx = 0; nodeIdx < numFineNodes; ++nodeIdx) {
316 if(PCrs->getNumEntriesInLocalRow(nodeIdx*dofsPerNode) == 1) {
319 for(LO dof = 0; dof < dofsPerNode; ++dof) {
320 rowIdx = nodeIdx*dofsPerNode + dof;
321 prolongatorGraph->getLocalRowView(rowIdx, colIndices);
322 PCrs->replaceLocalValues(rowIdx, colIndices, values());
328 for(
int dim = 0; dim < 3; ++dim) {
329 coords[0][dim] = fineCoords[dim][nodeIdx];
331 prolongatorGraph->getLocalRowView(nodeIdx*dofsPerNode, colIndices);
332 for(
int interpolationIdx=0; interpolationIdx < numInterpolationPoints; ++interpolationIdx) {
333 coords[interpolationIdx + 1].resize(3);
334 interpolationNodeIdx = colIndices[interpolationIdx] / dofsPerNode;
335 for(
int dim = 0; dim < 3; ++dim) {
336 coords[interpolationIdx + 1][dim] = ghostCoords[dim][interpolationNodeIdx];
339 ComputeLinearInterpolationStencil(numDimensions, numInterpolationPoints, coords, stencil);
340 values.resize(numInterpolationPoints);
341 for(LO valueIdx = 0; valueIdx < numInterpolationPoints; ++valueIdx) {
342 values[valueIdx] = Teuchos::as<SC>(stencil[valueIdx]);
346 for(LO dof = 0; dof < dofsPerNode; ++dof) {
347 rowIdx = nodeIdx*dofsPerNode + dof;
348 prolongatorGraph->getLocalRowView(rowIdx, colIndices);
349 PCrs->replaceLocalValues(rowIdx, colIndices, values());
354 *out <<
"The calculation of the interpolation stencils has completed." << std::endl;
356 PCrs->fillComplete();
358 *out <<
"All values in P have been set and expertStaticFillComplete has been performed." << std::endl;
361 if (A->IsView(
"stridedMaps") ==
true) {
362 P->CreateView(
"stridedMaps", A->getRowMap(
"stridedMaps"), stridedDomainMap);
364 P->CreateView(
"stridedMaps", P->getRangeMap(), stridedDomainMap);
370 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
373 const Array<Array<real_type> > coord,
374 Array<real_type>& stencil)
const {
391 Teuchos::SerialDenseMatrix<LO,real_type> Jacobian(numDimensions, numDimensions);
392 Teuchos::SerialDenseVector<LO,real_type> residual(numDimensions);
393 Teuchos::SerialDenseVector<LO,real_type> solutionDirection(numDimensions);
394 Teuchos::SerialDenseVector<LO,real_type> paramCoords(numDimensions);
395 Teuchos::SerialDenseSolver<LO,real_type> problem;
396 int iter = 0, max_iter = 5;
397 real_type functions[4][8], norm_ref = 1.0, norm2 = 1.0, tol = 1.0e-5;
398 paramCoords.size(numDimensions);
400 while( (iter < max_iter) && (norm2 > tol*norm_ref) ) {
403 solutionDirection.size(numDimensions);
404 residual.size(numDimensions);
408 GetInterpolationFunctions(numDimensions, paramCoords, functions);
409 for(LO i = 0; i < numDimensions; ++i) {
410 residual(i) = coord[0][i];
411 for(LO k = 0; k < numInterpolationPoints; ++k) {
412 residual(i) -= functions[0][k]*coord[k+1][i];
415 norm_ref += residual(i)*residual(i);
416 if(i == numDimensions - 1) {
417 norm_ref = std::sqrt(norm_ref);
421 for(LO j = 0; j < numDimensions; ++j) {
422 for(LO k = 0; k < numInterpolationPoints; ++k) {
423 Jacobian(i,j) += functions[j+1][k]*coord[k+1][i];
429 problem.setMatrix(Teuchos::rcp(&Jacobian,
false));
430 problem.setVectors(Teuchos::rcp(&solutionDirection,
false), Teuchos::rcp(&residual,
false));
431 if(problem.shouldEquilibrate()) {problem.factorWithEquilibration(
true);}
434 for(LO i = 0; i < numDimensions; ++i) {
435 paramCoords(i) = paramCoords(i) + solutionDirection(i);
439 GetInterpolationFunctions(numDimensions, paramCoords, functions);
440 for(LO i = 0; i < numDimensions; ++i) {
442 for(LO k = 0; k < numInterpolationPoints; ++k) {
443 tmp -= functions[0][k]*coord[k+1][i];
448 norm2 = std::sqrt(norm2);
452 for(LO i = 0; i < numInterpolationPoints; ++i) {
453 stencil[i] = functions[0][i];
458 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
461 const Teuchos::SerialDenseVector<LO, real_type> parametricCoordinates,
463 real_type xi = 0.0, eta = 0.0, zeta = 0.0, denominator = 0.0;
464 if(numDimensions == 1) {
465 xi = parametricCoordinates[0];
467 }
else if(numDimensions == 2) {
468 xi = parametricCoordinates[0];
469 eta = parametricCoordinates[1];
471 }
else if(numDimensions == 3) {
472 xi = parametricCoordinates[0];
473 eta = parametricCoordinates[1];
474 zeta = parametricCoordinates[2];
478 functions[0][0] = (1.0 - xi)*(1.0 - eta)*(1.0 - zeta) / denominator;
479 functions[0][1] = (1.0 + xi)*(1.0 - eta)*(1.0 - zeta) / denominator;
480 functions[0][2] = (1.0 - xi)*(1.0 + eta)*(1.0 - zeta) / denominator;
481 functions[0][3] = (1.0 + xi)*(1.0 + eta)*(1.0 - zeta) / denominator;
482 functions[0][4] = (1.0 - xi)*(1.0 - eta)*(1.0 + zeta) / denominator;
483 functions[0][5] = (1.0 + xi)*(1.0 - eta)*(1.0 + zeta) / denominator;
484 functions[0][6] = (1.0 - xi)*(1.0 + eta)*(1.0 + zeta) / denominator;
485 functions[0][7] = (1.0 + xi)*(1.0 + eta)*(1.0 + zeta) / denominator;
487 functions[1][0] = -(1.0 - eta)*(1.0 - zeta) / denominator;
488 functions[1][1] = (1.0 - eta)*(1.0 - zeta) / denominator;
489 functions[1][2] = -(1.0 + eta)*(1.0 - zeta) / denominator;
490 functions[1][3] = (1.0 + eta)*(1.0 - zeta) / denominator;
491 functions[1][4] = -(1.0 - eta)*(1.0 + zeta) / denominator;
492 functions[1][5] = (1.0 - eta)*(1.0 + zeta) / denominator;
493 functions[1][6] = -(1.0 + eta)*(1.0 + zeta) / denominator;
494 functions[1][7] = (1.0 + eta)*(1.0 + zeta) / denominator;
496 functions[2][0] = -(1.0 - xi)*(1.0 - zeta) / denominator;
497 functions[2][1] = -(1.0 + xi)*(1.0 - zeta) / denominator;
498 functions[2][2] = (1.0 - xi)*(1.0 - zeta) / denominator;
499 functions[2][3] = (1.0 + xi)*(1.0 - zeta) / denominator;
500 functions[2][4] = -(1.0 - xi)*(1.0 + zeta) / denominator;
501 functions[2][5] = -(1.0 + xi)*(1.0 + zeta) / denominator;
502 functions[2][6] = (1.0 - xi)*(1.0 + zeta) / denominator;
503 functions[2][7] = (1.0 + xi)*(1.0 + zeta) / denominator;
505 functions[3][0] = -(1.0 - xi)*(1.0 - eta) / denominator;
506 functions[3][1] = -(1.0 + xi)*(1.0 - eta) / denominator;
507 functions[3][2] = -(1.0 - xi)*(1.0 + eta) / denominator;
508 functions[3][3] = -(1.0 + xi)*(1.0 + eta) / denominator;
509 functions[3][4] = (1.0 - xi)*(1.0 - eta) / denominator;
510 functions[3][5] = (1.0 + xi)*(1.0 - eta) / denominator;
511 functions[3][6] = (1.0 - xi)*(1.0 + eta) / denominator;
512 functions[3][7] = (1.0 + xi)*(1.0 + eta) / denominator;
516 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
521 : geoData_(*geoData), fineCoordView_(fineCoordView), coarseCoordView_(coarseCoordView) {
525 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
526 KOKKOS_INLINE_FUNCTION
531 LO nodeCoarseTuple[3] = {0, 0, 0};
532 LO nodeFineTuple[3] = {0, 0, 0};
533 auto coarseningRate = geoData_.getCoarseningRates();
534 auto fineNodesPerDir = geoData_.getLocalFineNodesPerDir();
535 auto coarseNodesPerDir = geoData_.getCoarseNodesPerDir();
536 geoData_.getCoarseLID2CoarseTuple(coarseNodeIdx, nodeCoarseTuple);
537 for(
int dim = 0; dim < 3; ++dim) {
538 if(nodeCoarseTuple[dim] == coarseNodesPerDir(dim) - 1) {
539 nodeFineTuple[dim] = fineNodesPerDir(dim) - 1;
541 nodeFineTuple[dim] = nodeCoarseTuple[dim]*coarseningRate(dim);
545 fineNodeIdx = nodeFineTuple[2]*fineNodesPerDir(1)*fineNodesPerDir(0)
546 + nodeFineTuple[1]*fineNodesPerDir(0) + nodeFineTuple[0];
548 for(
int dim = 0; dim < fineCoordView_.extent_int(1); ++dim) {
549 coarseCoordView_(coarseNodeIdx, dim) = fineCoordView_(fineNodeIdx, dim);
555 #endif // MUELU_GEOMETRICINTERPOLATIONPFACTORY_KOKKOS_DEF_HPP
void GetInterpolationFunctions(const LO numDimensions, const Teuchos::SerialDenseVector< LO, real_type > parametricCoordinates, real_type functions[4][8]) const
coarseCoordinatesBuilderFunctor(RCP< IndexManager_kokkos > geoData, coord_view_type fineCoordView, coord_view_type coarseCoordView)
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
Timer to be used in factories. Similar to Monitor but with additional timers.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
Namespace for MueLu classes and methods.
typename Teuchos::ScalarTraits< SC >::coordinateType real_type
void BuildLinearP(RCP< Matrix > &A, RCP< const CrsGraph > &prolongatorGraph, RCP< realvaluedmultivector_type > &fineCoordinates, RCP< realvaluedmultivector_type > &ghostCoordinates, const int numDimensions, RCP< Matrix > &P) const
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Class that holds all level-specific information.
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
typename Kokkos::View< impl_scalar_type **, Kokkos::LayoutLeft, device_type > coord_view_type
#define SET_VALID_ENTRY(name)
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
void BuildConstantP(RCP< Matrix > &P, RCP< const CrsGraph > &prolongatorGraph, RCP< Matrix > &A) const
void ComputeLinearInterpolationStencil(const int numDimensions, const int numInterpolationPoints, const Array< Array< real_type > > coord, Array< real_type > &stencil) const
KOKKOS_INLINE_FUNCTION void operator()(const LO nodeIdx) const