75 const std::unordered_map<std::string,double>* elementBlockMultipliers)
77 PANZER_FUNC_TIME_MONITOR_DIFF(
"L2Projection::Build Mass Matrix",TopBuildMassMatrix);
78 TEUCHOS_ASSERT(setupCalled_);
80 if (elementBlockMultipliers !=
nullptr) {
81 TEUCHOS_ASSERT(elementBlockMultipliers->size() == elementBlockNames_.size());
85 std::vector<Teuchos::RCP<const panzer::GlobalIndexer>> indexers;
86 indexers.push_back(targetGlobalIndexer_);
92 ownedMatrix->resumeFill();
93 ownedMatrix->setAllToScalar(0.0);
94 ghostedMatrix->resumeFill();
95 ghostedMatrix->setAllToScalar(0.0);
97 const int fieldIndex = targetGlobalIndexer_->getFieldNum(targetBasisDescriptor_.getType());
99 const bool is_scalar = targetBasisDescriptor_.getType()==
"HGrad" || targetBasisDescriptor_.getType()==
"Const" || targetBasisDescriptor_.getType()==
"HVol";
103 auto M = ghostedMatrix->getLocalMatrixDevice();
104 for (
const auto& block : elementBlockNames_) {
106 double ebMultiplier = 1.0;
107 if (elementBlockMultipliers !=
nullptr)
108 ebMultiplier = elementBlockMultipliers->find(block)->second;
112 const auto worksets = worksetContainer_->getWorksets(wd);
114 for (
const auto& workset : *worksets) {
116 const auto basisValues = workset.getBasisValues(targetBasisDescriptor_,integrationDescriptor_);
118 const auto unweightedBasis = basisValues.getBasisValues(
false).get_static_view();
119 const auto weightedBasis = basisValues.getBasisValues(
true).get_static_view();
121 const std::vector<panzer::LocalOrdinal>&
offsets = targetGlobalIndexer_->getGIDFieldOffsets(block,fieldIndex);
122 PHX::View<panzer::LocalOrdinal*> kOffsets(
"MassMatrix: Offsets",
offsets.size());
123 auto kOffsets_h = Kokkos::create_mirror_view(kOffsets);
128 Kokkos::deep_copy(kOffsets, kOffsets_h);
131 PHX::View<panzer::LocalOrdinal**> localIds(
"MassMatrix: LocalIds", workset.numOwnedCells()+workset.numGhostCells()+workset.numVirtualCells(),
132 targetGlobalIndexer_->getElementBlockGIDCount(block));
135 const auto cellLocalIdsNoGhost = Kokkos::subview(workset.getLocalCellIDs(),std::make_pair(0,workset.numOwnedCells()));
137 targetGlobalIndexer_->getElementLIDs(cellLocalIdsNoGhost,localIds);
139 const int numBasisPoints =
static_cast<int>(weightedBasis.extent(1));
141 Kokkos::parallel_for(workset.numOwnedCells(),KOKKOS_LAMBDA (
const int& cell) {
142 double total_mass = 0.0, trace = 0.0;
144 panzer::LocalOrdinal cLIDs[256];
145 const int numIds =
static_cast<int>(localIds.extent(1));
146 for(
int i=0;i<numIds;++i)
147 cLIDs[i] = localIds(cell,i);
149 double vals[256]={0.0};
150 const int numQP =
static_cast<int>(unweightedBasis.extent(2));
152 for (
int row=0; row < numBasisPoints; ++row) {
153 for (
int col=0; col < numIds; ++col) {
154 for (
int qp=0; qp < numQP; ++qp) {
155 auto tmp = unweightedBasis(cell,row,qp) * weightedBasis(cell,col,qp) * ebMultiplier;
163 for (
int row=0; row < numBasisPoints; ++row) {
164 for (
int col=0; col < numBasisPoints; ++col)
167 int offset = kOffsets(row);
168 panzer::LocalOrdinal lid = localIds(cell,offset);
171 for (
int qp=0; qp < numQP; ++qp)
172 vals[col] += unweightedBasis(cell,row,qp) * weightedBasis(cell,col,qp) * ebMultiplier * total_mass / trace;
174 M.sumIntoValues(lid,cLIDs,numIds,
vals,
true,
true);
179 Kokkos::parallel_for(workset.numOwnedCells(),KOKKOS_LAMBDA (
const int& cell) {
181 panzer::LocalOrdinal cLIDs[256];
182 const int numIds =
static_cast<int>(localIds.extent(1));
183 for(
int i=0;i<numIds;++i)
184 cLIDs[i] = localIds(cell,i);
187 const int numQP =
static_cast<int>(unweightedBasis.extent(2));
189 for (
int row=0; row < numBasisPoints; ++row) {
190 int offset = kOffsets(row);
191 panzer::LocalOrdinal lid = localIds(cell,offset);
193 for (
int col=0; col < numIds; ++col) {
195 for (
int qp=0; qp < numQP; ++qp)
196 vals[col] += unweightedBasis(cell,row,qp) * weightedBasis(cell,col,qp) * ebMultiplier;
198 M.sumIntoValues(lid,cLIDs,numIds,
vals,
true,
true);
207 auto M = ghostedMatrix->getLocalMatrixDevice();
208 for (
const auto& block : elementBlockNames_) {
210 double ebMultiplier = 1.0;
211 if (elementBlockMultipliers !=
nullptr)
212 ebMultiplier = elementBlockMultipliers->find(block)->second;
216 const auto& worksets = worksetContainer_->getWorksets(wd);
218 for (
const auto& workset : *worksets) {
220 const auto basisValues = workset.getBasisValues(targetBasisDescriptor_,integrationDescriptor_);
222 const auto unweightedBasis = basisValues.getVectorBasisValues(
false).get_static_view();
223 const auto weightedBasis = basisValues.getVectorBasisValues(
true).get_static_view();
225 const std::vector<panzer::LocalOrdinal>&
offsets = targetGlobalIndexer_->getGIDFieldOffsets(block,fieldIndex);
226 PHX::View<panzer::LocalOrdinal*> kOffsets(
"MassMatrix: Offsets",
offsets.size());
227 auto kOffsets_h = Kokkos::create_mirror_view(kOffsets);
232 Kokkos::deep_copy(kOffsets, kOffsets_h);
235 PHX::View<panzer::LocalOrdinal**> localIds(
"MassMatrix: LocalIds", workset.numOwnedCells()+workset.numGhostCells()+workset.numVirtualCells(),
236 targetGlobalIndexer_->getElementBlockGIDCount(block));
239 const PHX::View<const int*> cellLocalIdsNoGhost = Kokkos::subview(workset.cell_local_ids_k,std::make_pair(0,workset.numOwnedCells()));
241 targetGlobalIndexer_->getElementLIDs(cellLocalIdsNoGhost,localIds);
243 const int numBasisPoints =
static_cast<int>(weightedBasis.extent(1));
244 Kokkos::parallel_for(workset.numOwnedCells(),KOKKOS_LAMBDA (
const int& cell) {
246 panzer::LocalOrdinal cLIDs[256];
247 const int numIds =
static_cast<int>(localIds.extent(1));
248 for(
int i=0;i<numIds;++i)
249 cLIDs[i] = localIds(cell,i);
252 const int numQP =
static_cast<int>(unweightedBasis.extent(2));
254 for (
int qp=0; qp < numQP; ++qp) {
255 for (
int row=0; row < numBasisPoints; ++row) {
256 int offset = kOffsets(row);
257 panzer::LocalOrdinal lid = localIds(cell,offset);
259 for (
int col=0; col < numIds; ++col){
261 for(
int dim=0; dim < static_cast<int>(weightedBasis.extent(3)); ++dim)
262 vals[col] += unweightedBasis(cell,row,qp,dim) * weightedBasis(cell,col,qp,dim) * ebMultiplier;
265 M.sumIntoValues(lid,cLIDs,numIds,
vals,
true,
true);
275 PANZER_FUNC_TIME_MONITOR_DIFF(
"Exporting of mass matrix",ExportMM);
276 auto map = factory.
getMap(0);
277 ghostedMatrix->fillComplete(map,map);
279 ownedMatrix->doExport(*ghostedMatrix, *exporter, Tpetra::ADD);
280 ownedMatrix->fillComplete();
307 const Teuchos::RCP<
const Tpetra::Map<panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType>>& inputOwnedSourceMap,
308 const std::string& sourceFieldName,
310 const int directionIndex)
317 using MapType = Tpetra::Map<panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType>;
318 using GraphType = Tpetra::CrsGraph<panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType>;
319 using ExportType = Tpetra::Export<panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType>;
320 using MatrixType = Tpetra::CrsMatrix<double,panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType>;
326 RCP<MapType> ghostedTargetMap;
328 std::vector<panzer::GlobalOrdinal> indices;
329 targetGlobalIndexer_->getOwnedAndGhostedIndices(indices);
330 ghostedTargetMap = rcp(
new MapType(Teuchos::OrdinalTraits<panzer::GlobalOrdinal>::invalid(),indices,0,comm_));
333 RCP<MapType> ghostedSourceMap;
335 std::vector<panzer::GlobalOrdinal> indices;
337 ghostedSourceMap = rcp(
new MapType(Teuchos::OrdinalTraits<panzer::GlobalOrdinal>::invalid(),indices,0,comm_));
343 std::vector<size_t> nEntriesPerRow(ghostedTargetMap->getLocalNumElements(),0);
344 std::vector<std::string> elementBlockIds;
345 targetGlobalIndexer_->getElementBlockIds(elementBlockIds);
346 std::vector<std::string>::const_iterator blockItr;
347 for (blockItr=elementBlockIds.begin();blockItr!=elementBlockIds.end();++blockItr) {
348 std::string blockId = *blockItr;
349 const std::vector<panzer::LocalOrdinal> & elements = targetGlobalIndexer_->getElementBlock(blockId);
351 std::vector<panzer::GlobalOrdinal> row_gids;
352 std::vector<panzer::GlobalOrdinal> col_gids;
354 for(std::size_t elmt=0;elmt<elements.size();elmt++) {
355 targetGlobalIndexer_->getElementGIDs(elements[elmt],row_gids);
357 for(std::size_t row=0;row<row_gids.size();row++) {
358 panzer::LocalOrdinal lid =
359 ghostedTargetMap->getLocalElement(row_gids[row]);
360 nEntriesPerRow[lid] += col_gids.size();
365 Teuchos::ArrayView<const size_t> nEntriesPerRowView(nEntriesPerRow);
366 RCP<GraphType> ghostedGraph = rcp(
new GraphType(ghostedTargetMap,ghostedSourceMap,nEntriesPerRowView));
368 for (blockItr=elementBlockIds.begin();blockItr!=elementBlockIds.end();++blockItr) {
369 std::string blockId = *blockItr;
370 const std::vector<panzer::LocalOrdinal> & elements = targetGlobalIndexer_->getElementBlock(blockId);
372 std::vector<panzer::GlobalOrdinal> row_gids;
373 std::vector<panzer::GlobalOrdinal> col_gids;
375 for(std::size_t elmt=0;elmt<elements.size();elmt++) {
376 targetGlobalIndexer_->getElementGIDs(elements[elmt],row_gids);
378 for(std::size_t row=0;row<row_gids.size();row++)
379 ghostedGraph->insertGlobalIndices(row_gids[row],col_gids);
383 RCP<MapType> ownedTargetMap;
385 std::vector<panzer::GlobalOrdinal> indices;
386 targetGlobalIndexer_->getOwnedIndices(indices);
387 ownedTargetMap = rcp(
new MapType(Teuchos::OrdinalTraits<panzer::GlobalOrdinal>::invalid(),indices,0,comm_));
390 RCP<const MapType> ownedSourceMap = inputOwnedSourceMap;
391 if (is_null(ownedSourceMap)) {
392 std::vector<panzer::GlobalOrdinal> indices;
394 ownedSourceMap = rcp(
new MapType(Teuchos::OrdinalTraits<panzer::GlobalOrdinal>::invalid(),indices,0,comm_));
398 ghostedGraph->fillComplete(ownedSourceMap,ownedTargetMap);
404 RCP<GraphType> ownedGraph = rcp(
new GraphType(ownedTargetMap,0));
405 RCP<const ExportType> exporter = rcp(
new ExportType(ghostedTargetMap,ownedTargetMap));
406 ownedGraph->doExport(*ghostedGraph, *exporter, Tpetra::INSERT);
407 ownedGraph->fillComplete(ownedSourceMap,ownedTargetMap);
409 RCP<MatrixType> ghostedMatrix = rcp(
new MatrixType(ghostedGraph));
410 RCP<MatrixType> ownedMatrix = rcp(
new MatrixType(ownedGraph));
414 ghostedMatrix->setAllToScalar(0.0);
415 ownedMatrix->setAllToScalar(0.0);
420 for (
const auto& block : elementBlockNames_) {
423 const auto& worksets = worksetContainer_->getWorksets(wd);
424 for (
const auto& workset : *worksets) {
427 const auto& targetBasisValues = workset.getBasisValues(targetBasisDescriptor_,integrationDescriptor_);
428 const auto& targetWeightedBasis = targetBasisValues.getBasisValues(
true).get_static_view();
431 const auto& sourceBasisValues = workset.getBasisValues(sourceBasisDescriptor,integrationDescriptor_);
432 PHX::View<const double***> sourceUnweightedScalarBasis;
433 PHX::View<const double****> sourceUnweightedVectorBasis;
434 bool useRankThreeBasis =
false;
435 if ( (sourceBasisDescriptor.
getType() ==
"HGrad") || (sourceBasisDescriptor.
getType() ==
"Const") || (sourceBasisDescriptor.
getType() ==
"HVol") ) {
436 if (directionIndex == -1) {
437 sourceUnweightedScalarBasis = sourceBasisValues.getBasisValues(
false).get_static_view();
438 useRankThreeBasis =
true;
441 sourceUnweightedVectorBasis = sourceBasisValues.getGradBasisValues(
false).get_static_view();
445 sourceUnweightedVectorBasis = sourceBasisValues.getVectorBasisValues(
false).get_static_view();
449 PHX::View<panzer::LocalOrdinal**> targetLocalIds(
"buildRHSMatrix: targetLocalIds", workset.numOwnedCells(),
450 targetGlobalIndexer_->getElementBlockGIDCount(block));
451 PHX::View<panzer::LocalOrdinal**> sourceLocalIds(
"buildRHSMatrix: sourceLocalIds", workset.numOwnedCells(),
455 const auto cellLocalIdsNoGhost = Kokkos::subview(workset.cell_local_ids_k,std::make_pair(0,workset.numOwnedCells()));
456 targetGlobalIndexer_->getElementLIDs(cellLocalIdsNoGhost,targetLocalIds);
457 sourceGlobalIndexer.
getElementLIDs(cellLocalIdsNoGhost,sourceLocalIds);
461 PHX::View<panzer::LocalOrdinal*> targetFieldOffsets;
463 const auto fieldIndex = targetGlobalIndexer_->getFieldNum(targetBasisDescriptor_.getType());
464 const std::vector<panzer::LocalOrdinal>&
offsets = targetGlobalIndexer_->getGIDFieldOffsets(block,fieldIndex);
465 targetFieldOffsets = PHX::View<panzer::LocalOrdinal*>(
"L2Projection:buildRHS:targetFieldOffsets",
offsets.size());
466 const auto hostOffsets = Kokkos::create_mirror_view(targetFieldOffsets);
467 for(
size_t i=0; i <
offsets.size(); ++i)
469 Kokkos::deep_copy(targetFieldOffsets,hostOffsets);
472 PHX::View<panzer::LocalOrdinal*> sourceFieldOffsets;
474 const auto fieldIndex = sourceGlobalIndexer.
getFieldNum(sourceFieldName);
476 TEUCHOS_TEST_FOR_EXCEPTION(
offsets.size() == 0, std::runtime_error,
477 "ERROR: panzer::L2Projection::buildRHSMatrix() - The source field, \""
478 << sourceFieldName <<
"\", does not exist in element block \""
480 sourceFieldOffsets = PHX::View<panzer::LocalOrdinal*>(
"L2Projection:buildRHS:sourceFieldOffsets",
offsets.size());
481 const auto hostOffsets = Kokkos::create_mirror_view(sourceFieldOffsets);
482 for(
size_t i=0; i <
offsets.size(); ++i)
484 Kokkos::deep_copy(sourceFieldOffsets,hostOffsets);
487 const auto localMatrix = ghostedMatrix->getLocalMatrixDevice();
488 const int numRows =
static_cast<int>(targetWeightedBasis.extent(1));
491 if (useRankThreeBasis) {
492 tmpNumCols =
static_cast<int>(sourceUnweightedScalarBasis.extent(1));
493 tmpNumQP =
static_cast<int>(sourceUnweightedScalarBasis.extent(2));
496 tmpNumCols =
static_cast<int>(sourceUnweightedVectorBasis.extent(1));
497 tmpNumQP =
static_cast<int>(sourceUnweightedVectorBasis.extent(2));
499 const int numCols = tmpNumCols;
500 const int numQP = tmpNumQP;
502 Kokkos::parallel_for(Kokkos::RangePolicy<PHX::Device>(0,workset.numOwnedCells()),KOKKOS_LAMBDA (
const int& cell) {
503 panzer::LocalOrdinal cLIDs[256];
505 for (
int row = 0; row < numRows; ++row) {
506 const int rowOffset = targetFieldOffsets(row);
507 const int rowLID = targetLocalIds(cell,rowOffset);
508 for (
int col = 0; col < numCols; ++col)
511 for (
int col = 0; col < numCols; ++col) {
512 for (
int qp = 0; qp < numQP; ++qp) {
513 const int colOffset = sourceFieldOffsets(col);
514 const int colLID = sourceLocalIds(cell,colOffset);
516 if (useRankThreeBasis)
517 vals[col] += sourceUnweightedScalarBasis(cell,col,qp) * targetWeightedBasis(cell,row,qp);
519 vals[col] += sourceUnweightedVectorBasis(cell,col,qp,directionIndex) * targetWeightedBasis(cell,row,qp);
522 localMatrix.sumIntoValues(rowLID,cLIDs,numCols,
vals,
true,
true);
528 ghostedMatrix->fillComplete(ownedSourceMap,ownedTargetMap);
529 ownedMatrix->resumeFill();
530 ownedMatrix->doExport(*ghostedMatrix,*exporter,Tpetra::ADD);
531 ownedMatrix->fillComplete(ownedSourceMap,ownedTargetMap);