MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_GeometricInterpolationPFactory_kokkos_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_KOKKOS_DEF_HPP
47#define MUELU_GEOMETRICINTERPOLATIONPFACTORY_KOKKOS_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_IndexManager_kokkos.hpp"
55
56#ifdef HAVE_MUELU_TPETRA
57#include "Xpetra_TpetraCrsMatrix.hpp"
58#endif
59
60
61// Including this one last ensure that the short names of the above headers are defined properly
63
64namespace MueLu {
65
66 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
68 RCP<ParameterList> validParamList = rcp(new ParameterList());
69
70#define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
71 SET_VALID_ENTRY("interp: build coarse coordinates");
72#undef SET_VALID_ENTRY
73
74 // general variables needed in GeometricInterpolationPFactory_kokkos
75 validParamList->set<RCP<const FactoryBase> >("A", Teuchos::null,
76 "Generating factory of the matrix A");
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> >("Nullspace", Teuchos::null,
82 "Fine level nullspace used to construct the coarse level nullspace.");
83 validParamList->set<RCP<const FactoryBase> >("numDimensions", Teuchos::null,
84 "Number of spacial dimensions in the problem.");
85 validParamList->set<RCP<const FactoryBase> >("lCoarseNodesPerDim", Teuchos::null,
86 "Number of nodes per spatial dimension on the coarse grid.");
87 validParamList->set<RCP<const FactoryBase> >("indexManager", Teuchos::null,
88 "The index manager associated with the local mesh.");
89 validParamList->set<RCP<const FactoryBase> >("structuredInterpolationOrder", Teuchos::null,
90 "Interpolation order for constructing the prolongator.");
91
92 return validParamList;
93 }
94
95 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
97 DeclareInput(Level& fineLevel, Level& coarseLevel) const {
98 const ParameterList& pL = GetParameterList();
99
100 Input(fineLevel, "A");
101 Input(fineLevel, "Nullspace");
102 Input(fineLevel, "numDimensions");
103 Input(fineLevel, "prolongatorGraph");
104 Input(fineLevel, "lCoarseNodesPerDim");
105 Input(fineLevel, "structuredInterpolationOrder");
106
107 if( pL.get<bool>("interp: build coarse coordinates") ||
108 Get<int>(fineLevel, "structuredInterpolationOrder") == 1) {
109 Input(fineLevel, "Coordinates");
110 Input(fineLevel, "indexManager");
111 }
112
113 }
114
115 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
117 Build(Level& fineLevel, Level &coarseLevel) const {
118 return BuildP(fineLevel, coarseLevel);
119 }
120
121 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
123 BuildP(Level &fineLevel, Level &coarseLevel) const {
124 FactoryMonitor m(*this, "BuildP", coarseLevel);
125
126 // Set debug outputs based on environment variable
127 RCP<Teuchos::FancyOStream> out;
128 if(const char* dbg = std::getenv("MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEBUG")) {
129 out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
130 out->setShowAllFrontMatter(false).setShowProcRank(true);
131 } else {
132 out = Teuchos::getFancyOStream(rcp(new Teuchos::oblackholestream()));
133 }
134
135 *out << "Starting GeometricInterpolationPFactory_kokkos::BuildP." << std::endl;
136
137 // Get inputs from the parameter list
138 const ParameterList& pL = GetParameterList();
139 const bool buildCoarseCoordinates = pL.get<bool>("interp: build coarse coordinates");
140 const int interpolationOrder = Get<int>(fineLevel, "structuredInterpolationOrder");
141 const int numDimensions = Get<int>(fineLevel, "numDimensions");
142
143 // Declared main input/outputs to be retrieved and placed on the fine resp. coarse level
144 RCP<Matrix> A = Get<RCP<Matrix> >(fineLevel, "A");
145 RCP<const CrsGraph> prolongatorGraph = Get<RCP<CrsGraph> >(fineLevel, "prolongatorGraph");
146 RCP<realvaluedmultivector_type> fineCoordinates, coarseCoordinates;
147 RCP<Matrix> P;
148
149 // Check if we need to build coarse coordinates as they are used if we construct
150 // a linear interpolation prolongator
151 if(buildCoarseCoordinates || (interpolationOrder == 1)) {
152 SubFactoryMonitor sfm(*this, "BuildCoordinates", coarseLevel);
153
154 // Extract data from fine level
155 RCP<IndexManager_kokkos> geoData = Get<RCP<IndexManager_kokkos> >(fineLevel, "indexManager");
156 fineCoordinates = Get< RCP<realvaluedmultivector_type> >(fineLevel, "Coordinates");
157
158 // Build coarse coordinates map/multivector
159 RCP<const Map> coarseCoordsMap = MapFactory::Build(fineCoordinates->getMap()->lib(),
160 Teuchos::OrdinalTraits<GO>::invalid(),
161 geoData->getNumCoarseNodes(),
162 fineCoordinates->getMap()->getIndexBase(),
163 fineCoordinates->getMap()->getComm());
164 coarseCoordinates = Xpetra::MultiVectorFactory<real_type,LO,GO,Node>::
165 Build(coarseCoordsMap, fineCoordinates->getNumVectors());
166
167 // Construct and launch functor to fill coarse coordinates values
168 // function should take a const view really
169 coarseCoordinatesBuilderFunctor myCoarseCoordinatesBuilder(geoData,
170 fineCoordinates-> getDeviceLocalView(Xpetra::Access::ReadWrite),
171 coarseCoordinates->getDeviceLocalView(Xpetra::Access::OverwriteAll));
172 Kokkos::parallel_for("GeometricInterpolation: build coarse coordinates",
173 Kokkos::RangePolicy<execution_space>(0, geoData->getNumCoarseNodes()),
174 myCoarseCoordinatesBuilder);
175
176 Set(coarseLevel, "Coordinates", coarseCoordinates);
177 }
178
179 *out << "Fine and coarse coordinates have been loaded from the fine level and set on the coarse level." << std::endl;
180
181 if(interpolationOrder == 0) {
182 SubFactoryMonitor sfm(*this, "BuildConstantP", coarseLevel);
183 // Compute the prolongator using piece-wise constant interpolation
184 BuildConstantP(P, prolongatorGraph, A);
185 } else if(interpolationOrder == 1) {
186 // Compute the prolongator using piece-wise linear interpolation
187 // First get all the required coordinates to compute the local part of P
188 RCP<realvaluedmultivector_type> ghostCoordinates
189 = Xpetra::MultiVectorFactory<real_type,LO,GO,NO>::Build(prolongatorGraph->getColMap(),
190 fineCoordinates->getNumVectors());
191 RCP<const Import> ghostImporter = ImportFactory::Build(coarseCoordinates->getMap(),
192 prolongatorGraph->getColMap());
193 ghostCoordinates->doImport(*coarseCoordinates, *ghostImporter, Xpetra::INSERT);
194
195 SubFactoryMonitor sfm(*this, "BuildLinearP", coarseLevel);
196 BuildLinearP(A, prolongatorGraph, fineCoordinates, ghostCoordinates, numDimensions, P);
197 }
198
199 *out << "The prolongator matrix has been built." << std::endl;
200
201 {
202 SubFactoryMonitor sfm(*this, "BuildNullspace", coarseLevel);
203 // Build the coarse nullspace
204 RCP<MultiVector> fineNullspace = Get< RCP<MultiVector> > (fineLevel, "Nullspace");
205 RCP<MultiVector> coarseNullspace = MultiVectorFactory::Build(P->getDomainMap(),
206 fineNullspace->getNumVectors());
207 P->apply(*fineNullspace, *coarseNullspace, Teuchos::TRANS, Teuchos::ScalarTraits<SC>::one(),
208 Teuchos::ScalarTraits<SC>::zero());
209 Set(coarseLevel, "Nullspace", coarseNullspace);
210 }
211
212 *out << "The coarse nullspace is constructed and set on the coarse level." << std::endl;
213
214 Array<LO> lNodesPerDir = Get<Array<LO> >(fineLevel, "lCoarseNodesPerDim");
215 Set(coarseLevel, "numDimensions", numDimensions);
216 Set(coarseLevel, "lNodesPerDim", lNodesPerDir);
217 Set(coarseLevel, "P", P);
218
219 *out << "GeometricInterpolationPFactory_kokkos::BuildP has completed." << std::endl;
220
221 } // BuildP
222
223 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
225 BuildConstantP(RCP<Matrix>& P, RCP<const CrsGraph>& prolongatorGraph, RCP<Matrix>& A) const {
226
227 // Set debug outputs based on environment variable
228 RCP<Teuchos::FancyOStream> out;
229 if(const char* dbg = std::getenv("MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEBUG")) {
230 out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
231 out->setShowAllFrontMatter(false).setShowProcRank(true);
232 } else {
233 out = Teuchos::getFancyOStream(rcp(new Teuchos::oblackholestream()));
234 }
235
236 *out << "BuildConstantP" << std::endl;
237
238 std::vector<size_t> strideInfo(1);
239 strideInfo[0] = A->GetFixedBlockSize();
240 RCP<const StridedMap> stridedDomainMap =
241 StridedMapFactory::Build(prolongatorGraph->getDomainMap(), strideInfo);
242
243 *out << "Call prolongator constructor" << std::endl;
244 using helpers=Xpetra::Helpers<Scalar,LocalOrdinal,GlobalOrdinal,Node>;
245 if(helpers::isTpetraBlockCrs(A)) {
246#ifdef HAVE_MUELU_TPETRA
247 LO NSDim = A->GetStorageBlockSize();
248
249 // Build the exploded Map
250 // FIXME: Should look at doing this on device
251 RCP<const Map> BlockMap = prolongatorGraph->getDomainMap();
252 Teuchos::ArrayView<const GO> block_dofs = BlockMap->getLocalElementList();
253 Teuchos::Array<GO> point_dofs(block_dofs.size()*NSDim);
254 for(LO i=0, ct=0; i<block_dofs.size(); i++) {
255 for(LO j=0; j<NSDim; j++) {
256 point_dofs[ct] = block_dofs[i]*NSDim + j;
257 ct++;
258 }
259 }
260
261 RCP<const Map> PointMap = MapFactory::Build(BlockMap->lib(),
262 BlockMap->getGlobalNumElements() *NSDim,
263 point_dofs(),
264 BlockMap->getIndexBase(),
265 BlockMap->getComm());
266 strideInfo[0] = A->GetFixedBlockSize();
267 RCP<const StridedMap> stridedPointMap = StridedMapFactory::Build(PointMap, strideInfo);
268
269 RCP<Xpetra::CrsMatrix<SC,LO,GO,NO> > P_xpetra = Xpetra::CrsMatrixFactory<SC,LO,GO,NO>::BuildBlock(prolongatorGraph, PointMap, A->getRangeMap(),NSDim);
270 RCP<Xpetra::TpetraBlockCrsMatrix<SC,LO,GO,NO> > P_tpetra = rcp_dynamic_cast<Xpetra::TpetraBlockCrsMatrix<SC,LO,GO,NO> >(P_xpetra);
271 if(P_tpetra.is_null()) throw std::runtime_error("BuildConstantP: Matrix factory did not return a Tpetra::BlockCrsMatrix");
272 RCP<CrsMatrixWrap> P_wrap = rcp(new CrsMatrixWrap(P_xpetra));
273
274 const LO stride = strideInfo[0]*strideInfo[0];
275 const LO in_stride = strideInfo[0];
276 typename CrsMatrix::local_graph_type localGraph = prolongatorGraph->getLocalGraphDevice();
277 auto rowptr = localGraph.row_map;
278 auto indices = localGraph.entries;
279 auto values = P_tpetra->getTpetra_BlockCrsMatrix()->getValuesDeviceNonConst();
280
281 using ISC = typename Tpetra::BlockCrsMatrix<SC,LO,GO,NO>::impl_scalar_type;
282 ISC one = Teuchos::ScalarTraits<ISC>::one();
283
284 const Kokkos::TeamPolicy<execution_space> policy(prolongatorGraph->getLocalNumRows(), 1);
285
286 Kokkos::parallel_for("MueLu:GeoInterpFact::BuildConstantP::fill", policy,
287 KOKKOS_LAMBDA(const typename Kokkos::TeamPolicy<execution_space>::member_type &thread) {
288 auto row = thread.league_rank();
289 for(LO j = (LO)rowptr[row]; j<(LO) rowptr[row+1]; j++) {
290 LO block_offset = j*stride;
291 for(LO k=0; k<in_stride; k++)
292 values[block_offset + k*(in_stride+1) ] = one;
293 }
294 });
295
296 P = P_wrap;
297 if (A->IsView("stridedMaps") == true) {
298 P->CreateView("stridedMaps", A->getRowMap("stridedMaps"), stridedPointMap);
299 }
300 else {
301 P->CreateView("stridedMaps", P->getRangeMap(), PointMap);
302 }
303
304#else
305 throw std::runtime_error("GeometricInteroplationFactory::BuildConstantP(): BlockCrs requires Tpetra");
306#endif
307
308 }
309 else {
310 // Create the prolongator matrix and its associated objects
311 RCP<ParameterList> dummyList = rcp(new ParameterList());
312 P = rcp(new CrsMatrixWrap(prolongatorGraph, dummyList));
313 RCP<CrsMatrix> PCrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
314 PCrs->setAllToScalar(1.0);
315 PCrs->fillComplete();
316
317 // set StridingInformation of P
318 if (A->IsView("stridedMaps") == true) {
319 P->CreateView("stridedMaps", A->getRowMap("stridedMaps"), stridedDomainMap);
320 } else {
321 P->CreateView("stridedMaps", P->getRangeMap(), stridedDomainMap);
322 }
323 }
324
325 } // BuildConstantP
326
327 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
329 BuildLinearP(RCP<Matrix>& A, RCP<const CrsGraph>& prolongatorGraph,
330 RCP<realvaluedmultivector_type>& fineCoordinates,
331 RCP<realvaluedmultivector_type>& ghostCoordinates,
332 const int numDimensions, RCP<Matrix>& P) const {
333
334 // Set debug outputs based on environment variable
335 RCP<Teuchos::FancyOStream> out;
336 if(const char* dbg = std::getenv("MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEBUG")) {
337 out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
338 out->setShowAllFrontMatter(false).setShowProcRank(true);
339 } else {
340 out = Teuchos::getFancyOStream(rcp(new Teuchos::oblackholestream()));
341 }
342
343 *out << "Entering BuildLinearP" << std::endl;
344
345 // Extract coordinates for interpolation stencil calculations
346 const LO numFineNodes = fineCoordinates->getLocalLength();
347 const LO numGhostNodes = ghostCoordinates->getLocalLength();
348 Array<ArrayRCP<const real_type> > fineCoords(3);
349 Array<ArrayRCP<const real_type> > ghostCoords(3);
350 const real_type realZero = Teuchos::as<real_type>(0.0);
351 ArrayRCP<real_type> fineZero(numFineNodes, realZero);
352 ArrayRCP<real_type> ghostZero(numGhostNodes, realZero);
353 for(int dim = 0; dim < 3; ++dim) {
354 if(dim < numDimensions) {
355 fineCoords[dim] = fineCoordinates->getData(dim);
356 ghostCoords[dim] = ghostCoordinates->getData(dim);
357 } else {
358 fineCoords[dim] = fineZero;
359 ghostCoords[dim] = ghostZero;
360 }
361 }
362
363 *out << "Coordinates extracted from the multivectors!" << std::endl;
364
365 // Compute 2^numDimensions using bit logic to avoid round-off errors
366 const int numInterpolationPoints = 1 << numDimensions;
367 const int dofsPerNode = A->GetFixedBlockSize();
368
369 std::vector<size_t> strideInfo(1);
370 strideInfo[0] = dofsPerNode;
371 RCP<const StridedMap> stridedDomainMap =
372 StridedMapFactory::Build(prolongatorGraph->getDomainMap(), strideInfo);
373
374 *out << "The maps of P have been computed" << std::endl;
375
376 RCP<ParameterList> dummyList = rcp(new ParameterList());
377 P = rcp(new CrsMatrixWrap(prolongatorGraph, dummyList));
378 RCP<CrsMatrix> PCrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
379 PCrs->resumeFill(); // The Epetra matrix is considered filled at this point.
380
381 LO interpolationNodeIdx = 0, rowIdx = 0;
382 ArrayView<const LO> colIndices;
383 Array<SC> values;
384 Array<Array<real_type> > coords(numInterpolationPoints + 1);
385 Array<real_type> stencil(numInterpolationPoints);
386 for(LO nodeIdx = 0; nodeIdx < numFineNodes; ++nodeIdx) {
387 if(PCrs->getNumEntriesInLocalRow(nodeIdx*dofsPerNode) == 1) {
388 values.resize(1);
389 values[0] = 1.0;
390 for(LO dof = 0; dof < dofsPerNode; ++dof) {
391 rowIdx = nodeIdx*dofsPerNode + dof;
392 prolongatorGraph->getLocalRowView(rowIdx, colIndices);
393 PCrs->replaceLocalValues(rowIdx, colIndices, values());
394 }
395 } else {
396 // Extract the coordinates associated with the current node
397 // and the neighboring coarse nodes
398 coords[0].resize(3);
399 for(int dim = 0; dim < 3; ++dim) {
400 coords[0][dim] = fineCoords[dim][nodeIdx];
401 }
402 prolongatorGraph->getLocalRowView(nodeIdx*dofsPerNode, colIndices);
403 for(int interpolationIdx=0; interpolationIdx < numInterpolationPoints; ++interpolationIdx) {
404 coords[interpolationIdx + 1].resize(3);
405 interpolationNodeIdx = colIndices[interpolationIdx] / dofsPerNode;
406 for(int dim = 0; dim < 3; ++dim) {
407 coords[interpolationIdx + 1][dim] = ghostCoords[dim][interpolationNodeIdx];
408 }
409 }
410 ComputeLinearInterpolationStencil(numDimensions, numInterpolationPoints, coords, stencil);
411 values.resize(numInterpolationPoints);
412 for(LO valueIdx = 0; valueIdx < numInterpolationPoints; ++valueIdx) {
413 values[valueIdx] = Teuchos::as<SC>(stencil[valueIdx]);
414 }
415
416 // Set values in all the rows corresponding to nodeIdx
417 for(LO dof = 0; dof < dofsPerNode; ++dof) {
418 rowIdx = nodeIdx*dofsPerNode + dof;
419 prolongatorGraph->getLocalRowView(rowIdx, colIndices);
420 PCrs->replaceLocalValues(rowIdx, colIndices, values());
421 }
422 }
423 }
424
425 *out << "The calculation of the interpolation stencils has completed." << std::endl;
426
427 PCrs->fillComplete();
428
429 *out << "All values in P have been set and expertStaticFillComplete has been performed." << std::endl;
430
431 // set StridingInformation of P
432 if (A->IsView("stridedMaps") == true) {
433 P->CreateView("stridedMaps", A->getRowMap("stridedMaps"), stridedDomainMap);
434 } else {
435 P->CreateView("stridedMaps", P->getRangeMap(), stridedDomainMap);
436 }
437
438 } // BuildLinearP
439
440
441 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
443 ComputeLinearInterpolationStencil(const int numDimensions, const int numInterpolationPoints,
444 const Array<Array<real_type> > coord,
445 Array<real_type>& stencil) const {
446
447 // 7 8 Find xi, eta and zeta such that
448 // x---------x
449 // /| /| Rx = x_p - sum N_i(xi,eta,zeta)x_i = 0
450 // 5/ | 6/ | Ry = y_p - sum N_i(xi,eta,zeta)y_i = 0
451 // x---------x | Rz = z_p - sum N_i(xi,eta,zeta)z_i = 0
452 // | | *P | |
453 // | x------|--x We can do this with a Newton solver:
454 // | /3 | /4 We will start with initial guess (xi,eta,zeta) = (0,0,0)
455 // |/ |/ We compute the Jacobian and iterate until convergence...
456 // z y x---------x
457 // | / 1 2 Once we have (xi,eta,zeta), we can evaluate all N_i which
458 // |/ give us the weights for the interpolation stencil!
459 // o---x
460 //
461
462 Teuchos::SerialDenseMatrix<LO,real_type> Jacobian(numDimensions, numDimensions);
463 Teuchos::SerialDenseVector<LO,real_type> residual(numDimensions);
464 Teuchos::SerialDenseVector<LO,real_type> solutionDirection(numDimensions);
465 Teuchos::SerialDenseVector<LO,real_type> paramCoords(numDimensions);
466 Teuchos::SerialDenseSolver<LO,real_type> problem;
467 int iter = 0, max_iter = 5;
468 real_type functions[4][8], norm_ref = 1.0, norm2 = 1.0, tol = 1.0e-5;
469 paramCoords.size(numDimensions);
470
471 while( (iter < max_iter) && (norm2 > tol*norm_ref) ) {
472 ++iter;
473 norm2 = 0.0;
474 solutionDirection.size(numDimensions);
475 residual.size(numDimensions);
476 Jacobian = 0.0;
477
478 // Compute Jacobian and Residual
479 GetInterpolationFunctions(numDimensions, paramCoords, functions);
480 for(LO i = 0; i < numDimensions; ++i) {
481 residual(i) = coord[0][i]; // Add coordinates from point of interest
482 for(LO k = 0; k < numInterpolationPoints; ++k) {
483 residual(i) -= functions[0][k]*coord[k+1][i]; //Remove contribution from all coarse points
484 }
485 if(iter == 1) {
486 norm_ref += residual(i)*residual(i);
487 if(i == numDimensions - 1) {
488 norm_ref = std::sqrt(norm_ref);
489 }
490 }
491
492 for(LO j = 0; j < numDimensions; ++j) {
493 for(LO k = 0; k < numInterpolationPoints; ++k) {
494 Jacobian(i,j) += functions[j+1][k]*coord[k+1][i];
495 }
496 }
497 }
498
499 // Set Jacobian, Vectors and solve problem
500 problem.setMatrix(Teuchos::rcp(&Jacobian, false));
501 problem.setVectors(Teuchos::rcp(&solutionDirection, false), Teuchos::rcp(&residual, false));
502 if(problem.shouldEquilibrate()) {problem.factorWithEquilibration(true);}
503 problem.solve();
504
505 for(LO i = 0; i < numDimensions; ++i) {
506 paramCoords(i) = paramCoords(i) + solutionDirection(i);
507 }
508
509 // Recompute Residual norm
510 GetInterpolationFunctions(numDimensions, paramCoords, functions);
511 for(LO i = 0; i < numDimensions; ++i) {
512 real_type tmp = coord[0][i];
513 for(LO k = 0; k < numInterpolationPoints; ++k) {
514 tmp -= functions[0][k]*coord[k+1][i];
515 }
516 norm2 += tmp*tmp;
517 tmp = 0.0;
518 }
519 norm2 = std::sqrt(norm2);
520 }
521
522 // Load the interpolation values onto the stencil.
523 for(LO i = 0; i < numInterpolationPoints; ++i) {
524 stencil[i] = functions[0][i];
525 }
526
527 } // End ComputeLinearInterpolationStencil
528
529 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
531 GetInterpolationFunctions(const LO numDimensions,
532 const Teuchos::SerialDenseVector<LO, real_type> parametricCoordinates,
533 real_type functions[4][8]) const {
534 real_type xi = 0.0, eta = 0.0, zeta = 0.0, denominator = 0.0;
535 if(numDimensions == 1) {
536 xi = parametricCoordinates[0];
537 denominator = 2.0;
538 } else if(numDimensions == 2) {
539 xi = parametricCoordinates[0];
540 eta = parametricCoordinates[1];
541 denominator = 4.0;
542 } else if(numDimensions == 3) {
543 xi = parametricCoordinates[0];
544 eta = parametricCoordinates[1];
545 zeta = parametricCoordinates[2];
546 denominator = 8.0;
547 }
548
549 functions[0][0] = (1.0 - xi)*(1.0 - eta)*(1.0 - zeta) / denominator;
550 functions[0][1] = (1.0 + xi)*(1.0 - eta)*(1.0 - zeta) / denominator;
551 functions[0][2] = (1.0 - xi)*(1.0 + eta)*(1.0 - zeta) / denominator;
552 functions[0][3] = (1.0 + xi)*(1.0 + eta)*(1.0 - zeta) / denominator;
553 functions[0][4] = (1.0 - xi)*(1.0 - eta)*(1.0 + zeta) / denominator;
554 functions[0][5] = (1.0 + xi)*(1.0 - eta)*(1.0 + zeta) / denominator;
555 functions[0][6] = (1.0 - xi)*(1.0 + eta)*(1.0 + zeta) / denominator;
556 functions[0][7] = (1.0 + xi)*(1.0 + eta)*(1.0 + zeta) / denominator;
557
558 functions[1][0] = -(1.0 - eta)*(1.0 - zeta) / denominator;
559 functions[1][1] = (1.0 - eta)*(1.0 - zeta) / denominator;
560 functions[1][2] = -(1.0 + eta)*(1.0 - zeta) / denominator;
561 functions[1][3] = (1.0 + eta)*(1.0 - zeta) / denominator;
562 functions[1][4] = -(1.0 - eta)*(1.0 + zeta) / denominator;
563 functions[1][5] = (1.0 - eta)*(1.0 + zeta) / denominator;
564 functions[1][6] = -(1.0 + eta)*(1.0 + zeta) / denominator;
565 functions[1][7] = (1.0 + eta)*(1.0 + zeta) / denominator;
566
567 functions[2][0] = -(1.0 - xi)*(1.0 - zeta) / denominator;
568 functions[2][1] = -(1.0 + xi)*(1.0 - zeta) / denominator;
569 functions[2][2] = (1.0 - xi)*(1.0 - zeta) / denominator;
570 functions[2][3] = (1.0 + xi)*(1.0 - zeta) / denominator;
571 functions[2][4] = -(1.0 - xi)*(1.0 + zeta) / denominator;
572 functions[2][5] = -(1.0 + xi)*(1.0 + zeta) / denominator;
573 functions[2][6] = (1.0 - xi)*(1.0 + zeta) / denominator;
574 functions[2][7] = (1.0 + xi)*(1.0 + zeta) / denominator;
575
576 functions[3][0] = -(1.0 - xi)*(1.0 - eta) / denominator;
577 functions[3][1] = -(1.0 + xi)*(1.0 - eta) / denominator;
578 functions[3][2] = -(1.0 - xi)*(1.0 + eta) / denominator;
579 functions[3][3] = -(1.0 + xi)*(1.0 + eta) / denominator;
580 functions[3][4] = (1.0 - xi)*(1.0 - eta) / denominator;
581 functions[3][5] = (1.0 + xi)*(1.0 - eta) / denominator;
582 functions[3][6] = (1.0 - xi)*(1.0 + eta) / denominator;
583 functions[3][7] = (1.0 + xi)*(1.0 + eta) / denominator;
584
585 } // End GetInterpolationFunctions
586
587 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
590 coord_view_type fineCoordView,
591 coord_view_type coarseCoordView)
592 : geoData_(*geoData), fineCoordView_(fineCoordView), coarseCoordView_(coarseCoordView) {
593
594 } // coarseCoordinatesBuilderFunctor()
595
596 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
597 KOKKOS_INLINE_FUNCTION
599 coarseCoordinatesBuilderFunctor::operator() (const LO coarseNodeIdx) const {
600
601 LO fineNodeIdx;
602 LO nodeCoarseTuple[3] = {0, 0, 0};
603 LO nodeFineTuple[3] = {0, 0, 0};
604 auto coarseningRate = geoData_.getCoarseningRates();
605 auto fineNodesPerDir = geoData_.getLocalFineNodesPerDir();
606 auto coarseNodesPerDir = geoData_.getCoarseNodesPerDir();
607 geoData_.getCoarseLID2CoarseTuple(coarseNodeIdx, nodeCoarseTuple);
608 for(int dim = 0; dim < 3; ++dim) {
609 if(nodeCoarseTuple[dim] == coarseNodesPerDir(dim) - 1) {
610 nodeFineTuple[dim] = fineNodesPerDir(dim) - 1;
611 } else {
612 nodeFineTuple[dim] = nodeCoarseTuple[dim]*coarseningRate(dim);
613 }
614 }
615
616 fineNodeIdx = nodeFineTuple[2]*fineNodesPerDir(1)*fineNodesPerDir(0)
617 + nodeFineTuple[1]*fineNodesPerDir(0) + nodeFineTuple[0];
618
619 for(int dim = 0; dim < fineCoordView_.extent_int(1); ++dim) {
620 coarseCoordView_(coarseNodeIdx, dim) = fineCoordView_(fineNodeIdx, dim);
621 }
622 }
623
624} // namespace MueLu
625
626#endif // MUELU_GEOMETRICINTERPOLATIONPFACTORY_KOKKOS_DEF_HPP
#define SET_VALID_ENTRY(name)
Timer to be used in factories. Similar to Monitor but with additional timers.
typename Teuchos::ScalarTraits< SC >::coordinateType real_type
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
void ComputeLinearInterpolationStencil(const int numDimensions, const int numInterpolationPoints, const Array< Array< real_type > > coord, Array< real_type > &stencil) const
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.
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
typename Kokkos::View< impl_scalar_type **, Kokkos::LayoutLeft, device_type > coord_view_type
void BuildConstantP(RCP< Matrix > &P, RCP< const CrsGraph > &prolongatorGraph, RCP< Matrix > &A) const
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.
Definition: MueLu_Level.hpp:99
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
Namespace for MueLu classes and methods.
coarseCoordinatesBuilderFunctor(RCP< IndexManager_kokkos > geoData, coord_view_type fineCoordView, coord_view_type coarseCoordView)