46 #ifndef MUELU_GEOMETRICINTERPOLATIONPFACTORY_KOKKOS_DEF_HPP
47 #define MUELU_GEOMETRICINTERPOLATIONPFACTORY_KOKKOS_DEF_HPP
54 #include "MueLu_IndexManager_kokkos.hpp"
61 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
65 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
68 #undef SET_VALID_ENTRY
72 "Generating factory of the matrix A");
74 "Graph generated by StructuredAggregationFactory used to construct a piece-linear prolongator.");
76 "Fine level coordinates used to construct piece-wise linear prolongator and coarse level coordinates.");
78 "Fine level nullspace used to construct the coarse level nullspace.");
80 "Number of spacial dimensions in the problem.");
82 "Number of nodes per spatial dimension on the coarse grid.");
84 "The index manager associated with the local mesh.");
86 return validParamList;
89 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
94 Input(fineLevel,
"A");
95 Input(fineLevel,
"Nullspace");
96 Input(fineLevel,
"numDimensions");
97 Input(fineLevel,
"prolongatorGraph");
98 Input(fineLevel,
"lCoarseNodesPerDim");
100 if( pL.
get<
bool>(
"interp: build coarse coordinates") ||
101 (pL.
get<
int>(
"interp: interpolation order") == 1) ) {
102 Input(fineLevel,
"Coordinates");
103 Input(fineLevel,
"indexManager");
108 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
111 return BuildP(fineLevel, coarseLevel);
114 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
121 if(
const char* dbg = std::getenv(
"MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEBUG")) {
122 out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
128 *out <<
"Starting GeometricInterpolationPFactory_kokkos::BuildP." << std::endl;
132 const bool buildCoarseCoordinates = pL.
get<
bool>(
"interp: build coarse coordinates");
133 const int interpolationOrder = pL.
get<
int> (
"interp: interpolation order");
134 const int numDimensions = Get<int>(fineLevel,
"numDimensions");
144 if(buildCoarseCoordinates || (interpolationOrder == 1)) {
149 fineCoordinates = Get< RCP<realvaluedmultivector_type> >(fineLevel,
"Coordinates");
152 RCP<const Map> coarseCoordsMap = MapFactory::Build(fineCoordinates->getMap()->lib(),
154 geoData->getNumCoarseNodes(),
155 fineCoordinates->getMap()->getIndexBase(),
156 fineCoordinates->getMap()->getComm());
158 Build(coarseCoordsMap, fineCoordinates->getNumVectors());
162 fineCoordinates->
template getLocalView<device_type>(),
163 coarseCoordinates->template getLocalView<device_type>());
166 myCoarseCoordinatesBuilder);
168 Set(coarseLevel,
"Coordinates", coarseCoordinates);
171 *out <<
"Fine and coarse coordinates have been loaded from the fine level and set on the coarse level." << std::endl;
173 if(interpolationOrder == 0) {
176 BuildConstantP(P, prolongatorGraph, A);
177 }
else if(interpolationOrder == 1) {
182 fineCoordinates->getNumVectors());
183 RCP<const Import> ghostImporter = ImportFactory::Build(coarseCoordinates->getMap(),
184 prolongatorGraph->getColMap());
185 ghostCoordinates->doImport(*coarseCoordinates, *ghostImporter,
Xpetra::INSERT);
188 BuildLinearP(A, prolongatorGraph, fineCoordinates, ghostCoordinates, numDimensions, P);
191 *out <<
"The prolongator matrix has been built." << std::endl;
196 RCP<MultiVector> fineNullspace = Get< RCP<MultiVector> > (fineLevel,
"Nullspace");
197 RCP<MultiVector> coarseNullspace = MultiVectorFactory::Build(P->getDomainMap(),
198 fineNullspace->getNumVectors());
201 Set(coarseLevel,
"Nullspace", coarseNullspace);
204 *out <<
"The coarse nullspace is constructed and set on the coarse level." << std::endl;
206 Array<LO> lNodesPerDir = Get<Array<LO> >(fineLevel,
"lCoarseNodesPerDim");
207 Set(coarseLevel,
"numDimensions", numDimensions);
208 Set(coarseLevel,
"lNodesPerDim", lNodesPerDir);
209 Set(coarseLevel,
"P", P);
211 *out <<
"GeometricInterpolationPFactory_kokkos::BuildP has completed." << std::endl;
215 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
221 if(
const char* dbg = std::getenv(
"MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEBUG")) {
222 out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
228 *out <<
"BuildConstantP" << std::endl;
230 std::vector<size_t> strideInfo(1);
231 strideInfo[0] = A->GetFixedBlockSize();
233 StridedMapFactory::Build(prolongatorGraph->getDomainMap(), strideInfo);
235 *out <<
"Call prolongator constructor" << std::endl;
239 P =
rcp(
new CrsMatrixWrap(prolongatorGraph, dummyList));
240 RCP<CrsMatrix> PCrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
241 PCrs->setAllToScalar(1.0);
242 PCrs->fillComplete();
245 if (A->IsView(
"stridedMaps") ==
true) {
246 P->CreateView(
"stridedMaps", A->getRowMap(
"stridedMaps"), stridedDomainMap);
248 P->CreateView(
"stridedMaps", P->getRangeMap(), stridedDomainMap);
253 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
262 if(
const char* dbg = std::getenv(
"MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEBUG")) {
263 out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
269 *out <<
"Entering BuildLinearP" << std::endl;
272 const LO numFineNodes = fineCoordinates->getLocalLength();
273 const LO numGhostNodes = ghostCoordinates->getLocalLength();
276 const real_type realZero = Teuchos::as<real_type>(0.0);
279 for(
int dim = 0; dim < 3; ++dim) {
280 if(dim < numDimensions) {
281 fineCoords[dim] = fineCoordinates->getData(dim);
282 ghostCoords[dim] = ghostCoordinates->getData(dim);
284 fineCoords[dim] = fineZero;
285 ghostCoords[dim] = ghostZero;
289 *out <<
"Coordinates extracted from the multivectors!" << std::endl;
292 const int numInterpolationPoints = 1 << numDimensions;
293 const int dofsPerNode = A->GetFixedBlockSize();
295 std::vector<size_t> strideInfo(1);
296 strideInfo[0] = dofsPerNode;
298 StridedMapFactory::Build(prolongatorGraph->getDomainMap(), strideInfo);
300 *out <<
"The maps of P have been computed" << std::endl;
303 P =
rcp(
new CrsMatrixWrap(prolongatorGraph, dummyList));
304 RCP<CrsMatrix> PCrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
307 LO interpolationNodeIdx = 0, rowIdx = 0;
312 for(
LO nodeIdx = 0; nodeIdx < numFineNodes; ++nodeIdx) {
313 if(PCrs->getNumEntriesInLocalRow(nodeIdx*dofsPerNode) == 1) {
316 for(
LO dof = 0; dof < dofsPerNode; ++dof) {
317 rowIdx = nodeIdx*dofsPerNode + dof;
318 prolongatorGraph->getLocalRowView(rowIdx, colIndices);
319 PCrs->replaceLocalValues(rowIdx, colIndices, values());
325 for(
int dim = 0; dim < 3; ++dim) {
326 coords[0][dim] = fineCoords[dim][nodeIdx];
328 prolongatorGraph->getLocalRowView(nodeIdx*dofsPerNode, colIndices);
329 for(
int interpolationIdx=0; interpolationIdx < numInterpolationPoints; ++interpolationIdx) {
330 coords[interpolationIdx + 1].
resize(3);
331 interpolationNodeIdx = colIndices[interpolationIdx] / dofsPerNode;
332 for(
int dim = 0; dim < 3; ++dim) {
333 coords[interpolationIdx + 1][dim] = ghostCoords[dim][interpolationNodeIdx];
336 ComputeLinearInterpolationStencil(numDimensions, numInterpolationPoints, coords, stencil);
337 values.
resize(numInterpolationPoints);
338 for(
LO valueIdx = 0; valueIdx < numInterpolationPoints; ++valueIdx) {
339 values[valueIdx] = Teuchos::as<SC>(stencil[valueIdx]);
343 for(
LO dof = 0; dof < dofsPerNode; ++dof) {
344 rowIdx = nodeIdx*dofsPerNode + dof;
345 prolongatorGraph->getLocalRowView(rowIdx, colIndices);
346 PCrs->replaceLocalValues(rowIdx, colIndices, values());
351 *out <<
"The calculation of the interpolation stencils has completed." << std::endl;
353 PCrs->fillComplete();
355 *out <<
"All values in P have been set and expertStaticFillComplete has been performed." << std::endl;
358 if (A->IsView(
"stridedMaps") ==
true) {
359 P->CreateView(
"stridedMaps", A->getRowMap(
"stridedMaps"), stridedDomainMap);
361 P->CreateView(
"stridedMaps", P->getRangeMap(), stridedDomainMap);
367 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
393 int iter = 0, max_iter = 5;
394 real_type functions[4][8], norm_ref = 1.0, norm2 = 1.0, tol = 1.0e-5;
395 paramCoords.
size(numDimensions);
397 while( (iter < max_iter) && (norm2 > tol*norm_ref) ) {
400 solutionDirection.
size(numDimensions);
401 residual.
size(numDimensions);
405 GetInterpolationFunctions(numDimensions, paramCoords, functions);
406 for(
LO i = 0; i < numDimensions; ++i) {
407 residual(i) = coord[0][i];
408 for(
LO k = 0; k < numInterpolationPoints; ++k) {
409 residual(i) -= functions[0][k]*coord[k+1][i];
412 norm_ref += residual(i)*residual(i);
413 if(i == numDimensions - 1) {
414 norm_ref = std::sqrt(norm_ref);
418 for(
LO j = 0; j < numDimensions; ++j) {
419 for(
LO k = 0; k < numInterpolationPoints; ++k) {
420 Jacobian(i,j) += functions[j+1][k]*coord[k+1][i];
431 for(
LO i = 0; i < numDimensions; ++i) {
432 paramCoords(i) = paramCoords(i) + solutionDirection(i);
436 GetInterpolationFunctions(numDimensions, paramCoords, functions);
437 for(
LO i = 0; i < numDimensions; ++i) {
439 for(
LO k = 0; k < numInterpolationPoints; ++k) {
440 tmp -= functions[0][k]*coord[k+1][i];
445 norm2 = std::sqrt(norm2);
449 for(
LO i = 0; i < numInterpolationPoints; ++i) {
450 stencil[i] = functions[0][i];
455 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
460 real_type xi = 0.0, eta = 0.0, zeta = 0.0, denominator = 0.0;
461 if(numDimensions == 1) {
462 xi = parametricCoordinates[0];
464 }
else if(numDimensions == 2) {
465 xi = parametricCoordinates[0];
466 eta = parametricCoordinates[1];
468 }
else if(numDimensions == 3) {
469 xi = parametricCoordinates[0];
470 eta = parametricCoordinates[1];
471 zeta = parametricCoordinates[2];
475 functions[0][0] = (1.0 - xi)*(1.0 - eta)*(1.0 - zeta) / denominator;
476 functions[0][1] = (1.0 + xi)*(1.0 - eta)*(1.0 - zeta) / denominator;
477 functions[0][2] = (1.0 - xi)*(1.0 + eta)*(1.0 - zeta) / denominator;
478 functions[0][3] = (1.0 + xi)*(1.0 + eta)*(1.0 - zeta) / denominator;
479 functions[0][4] = (1.0 - xi)*(1.0 - eta)*(1.0 + zeta) / denominator;
480 functions[0][5] = (1.0 + xi)*(1.0 - eta)*(1.0 + zeta) / denominator;
481 functions[0][6] = (1.0 - xi)*(1.0 + eta)*(1.0 + zeta) / denominator;
482 functions[0][7] = (1.0 + xi)*(1.0 + eta)*(1.0 + zeta) / denominator;
484 functions[1][0] = -(1.0 - eta)*(1.0 - zeta) / denominator;
485 functions[1][1] = (1.0 - eta)*(1.0 - zeta) / denominator;
486 functions[1][2] = -(1.0 + eta)*(1.0 - zeta) / denominator;
487 functions[1][3] = (1.0 + eta)*(1.0 - zeta) / denominator;
488 functions[1][4] = -(1.0 - eta)*(1.0 + zeta) / denominator;
489 functions[1][5] = (1.0 - eta)*(1.0 + zeta) / denominator;
490 functions[1][6] = -(1.0 + eta)*(1.0 + zeta) / denominator;
491 functions[1][7] = (1.0 + eta)*(1.0 + zeta) / denominator;
493 functions[2][0] = -(1.0 - xi)*(1.0 - zeta) / denominator;
494 functions[2][1] = -(1.0 + xi)*(1.0 - zeta) / denominator;
495 functions[2][2] = (1.0 - xi)*(1.0 - zeta) / denominator;
496 functions[2][3] = (1.0 + xi)*(1.0 - zeta) / denominator;
497 functions[2][4] = -(1.0 - xi)*(1.0 + zeta) / denominator;
498 functions[2][5] = -(1.0 + xi)*(1.0 + zeta) / denominator;
499 functions[2][6] = (1.0 - xi)*(1.0 + zeta) / denominator;
500 functions[2][7] = (1.0 + xi)*(1.0 + zeta) / denominator;
502 functions[3][0] = -(1.0 - xi)*(1.0 - eta) / denominator;
503 functions[3][1] = -(1.0 + xi)*(1.0 - eta) / denominator;
504 functions[3][2] = -(1.0 - xi)*(1.0 + eta) / denominator;
505 functions[3][3] = -(1.0 + xi)*(1.0 + eta) / denominator;
506 functions[3][4] = (1.0 - xi)*(1.0 - eta) / denominator;
507 functions[3][5] = (1.0 + xi)*(1.0 - eta) / denominator;
508 functions[3][6] = (1.0 - xi)*(1.0 + eta) / denominator;
509 functions[3][7] = (1.0 + xi)*(1.0 + eta) / denominator;
513 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
518 : geoData_(*geoData), fineCoordView_(fineCoordView), coarseCoordView_(coarseCoordView) {
522 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
523 KOKKOS_INLINE_FUNCTION
528 LO nodeCoarseTuple[3] = {0, 0, 0};
529 LO nodeFineTuple[3] = {0, 0, 0};
530 auto coarseningRate = geoData_.getCoarseningRates();
531 auto fineNodesPerDir = geoData_.getLocalFineNodesPerDir();
532 auto coarseNodesPerDir = geoData_.getCoarseNodesPerDir();
533 geoData_.getCoarseLID2CoarseTuple(coarseNodeIdx, nodeCoarseTuple);
534 for(
int dim = 0; dim < 3; ++dim) {
535 if(nodeCoarseTuple[dim] == coarseNodesPerDir(dim) - 1) {
536 nodeFineTuple[dim] = fineNodesPerDir(dim) - 1;
538 nodeFineTuple[dim] = nodeCoarseTuple[dim]*coarseningRate(dim);
542 fineNodeIdx = nodeFineTuple[2]*fineNodesPerDir(1)*fineNodesPerDir(0)
543 + nodeFineTuple[1]*fineNodesPerDir(0) + nodeFineTuple[0];
545 for(
int dim = 0; dim < fineCoordView_.extent_int(1); ++dim) {
546 coarseCoordView_(coarseNodeIdx, dim) = fineCoordView_(fineNodeIdx, dim);
552 #endif // MUELU_GEOMETRICINTERPOLATIONPFACTORY_KOKKOS_DEF_HPP
void parallel_for(const ExecPolicy &policy, const FunctorType &functor, const std::string &str="", typename Impl::enable_if< Kokkos::Impl::is_execution_policy< ExecPolicy >::value >::type *=0)
basic_FancyOStream & setShowProcRank(const bool showProcRank)
coarseCoordinatesBuilderFunctor(RCP< IndexManager_kokkos > geoData, coord_view_type fineCoordView, coord_view_type coarseCoordView)
typename Kokkos::View< impl_scalar_type **, Kokkos::LayoutLeft, execution_space > coord_view_type
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Timer to be used in factories. Similar to Monitor but with additional timers.
void GetInterpolationFunctions(const LO numDimensions, const Teuchos::SerialDenseVector< LO, real_type > parametricCoordinates, real_type functions[4][8]) const
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
typename Teuchos::ScalarTraits< SC >::coordinateType real_type
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Class that holds all level-specific information.
void ComputeLinearInterpolationStencil(const int numDimensions, const int numInterpolationPoints, const Array< Array< real_type > > coord, Array< real_type > &stencil) const
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
void factorWithEquilibration(bool flag)
void BuildConstantP(RCP< Matrix > &P, RCP< const CrsGraph > &prolongatorGraph, RCP< Matrix > &A) const
#define SET_VALID_ENTRY(name)
basic_FancyOStream & setShowAllFrontMatter(const bool showAllFrontMatter)
void resize(size_type new_size, const value_type &x=value_type())
int size(OrdinalType length_in)
int setVectors(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &X, const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &B)
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
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.
static Teuchos::RCP< MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Build(const Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &map, size_t NumVectors, bool zeroOut=true)
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
KOKKOS_INLINE_FUNCTION void operator()(const LO nodeIdx) const
int setMatrix(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &A)