Panzer Version of the Day
Loading...
Searching...
No Matches
Panzer_L2Projection.cpp
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// Panzer: A partial differential equation assembly
4// engine for strongly coupled complex multiphysics systems
5//
6// Copyright 2011 NTESS and the Panzer contributors.
7// SPDX-License-Identifier: BSD-3-Clause
8// *****************************************************************************
9// @HEADER
10
11#ifndef PANZER_L2_PROJECTION_IMPL_HPP
12#define PANZER_L2_PROJECTION_IMPL_HPP
13
14#include "Teuchos_DefaultMpiComm.hpp"
15#include "Tpetra_CrsGraph.hpp"
16#include "Tpetra_MultiVector.hpp"
17#include "Tpetra_CrsMatrix.hpp"
18#include "Panzer_L2Projection.hpp"
19#include "Panzer_BasisDescriptor.hpp"
20#include "Panzer_TpetraLinearObjContainer.hpp"
21#include "Panzer_TpetraLinearObjFactory.hpp"
22#include "Panzer_BlockedTpetraLinearObjFactory.hpp"
23#include "Panzer_ElementBlockIdToPhysicsIdMap.hpp"
24#include "Panzer_DOFManagerFactory.hpp"
25#include "Panzer_BlockedDOFManagerFactory.hpp"
26#include "Panzer_IntrepidFieldPattern.hpp"
27#include "Panzer_WorksetDescriptor.hpp"
28#include "Panzer_WorksetContainer.hpp"
29#include "Panzer_Workset.hpp"
30
31namespace panzer {
32
33 void panzer::L2Projection::setup(const panzer::BasisDescriptor& targetBasis,
34 const panzer::IntegrationDescriptor& integrationDescriptor,
35 const Teuchos::RCP<const Teuchos::MpiComm<int>>& comm,
36 const Teuchos::RCP<const panzer::ConnManager>& connManager,
37 const std::vector<std::string>& elementBlockNames,
38 const Teuchos::RCP<panzer::WorksetContainer> worksetContainer)
39 {
40 targetBasisDescriptor_ = targetBasis;
41 integrationDescriptor_ = integrationDescriptor;
42 comm_ = comm;
43 connManager_ = connManager;
44 elementBlockNames_ = elementBlockNames;
45 worksetContainer_ = worksetContainer;
46 setupCalled_ = true;
47 useUserSuppliedBasisValues_ = false;
48
49 // Build target DOF Manager
50 targetGlobalIndexer_ =
51 Teuchos::rcp(new panzer::DOFManager(Teuchos::rcp_const_cast<panzer::ConnManager>(connManager),*(comm->getRawMpiComm())));
52
53 // For hybrid mesh, blocks could have different topology
54 for (const auto& eBlock : elementBlockNames_) {
55 std::vector<shards::CellTopology> topologies;
56 connManager_->getElementBlockTopologies(topologies);
57 std::vector<std::string> ebNames;
58 connManager_->getElementBlockIds(ebNames);
59 const auto search = std::find(ebNames.cbegin(),ebNames.cend(),eBlock);
60 TEUCHOS_ASSERT(search != ebNames.cend());
61 const int index = std::distance(ebNames.cbegin(),search);
62 const auto& cellTopology = topologies[index];
63
64 auto intrepidBasis = panzer::createIntrepid2Basis<PHX::Device,double,double>(targetBasisDescriptor_.getType(),
65 targetBasisDescriptor_.getOrder(),
66 cellTopology);
67 Teuchos::RCP<const panzer::FieldPattern> fieldPattern(new panzer::Intrepid2FieldPattern(intrepidBasis));
68 // Field name is the basis type
69 targetGlobalIndexer_->addField(eBlock,targetBasisDescriptor_.getType(),fieldPattern);
70 }
71
72 targetGlobalIndexer_->buildGlobalUnknowns();
73
74 // Check workset needs are correct
75 }
76
77 void panzer::L2Projection::useBasisValues(const std::map<std::string,Teuchos::RCP<panzer::BasisValues2<double>>>& bv)
78 {
79 useUserSuppliedBasisValues_ = true;
80 basisValues_ = bv;
81
82 // Check that the basis and integration descriptor match what was
83 // supplied in setup.
84 for (const auto& eblock : elementBlockNames_) {
85 TEUCHOS_ASSERT(basisValues_[eblock]->getBasisDescriptor()==targetBasisDescriptor_);
86 }
87 }
88
89 Teuchos::RCP<panzer::GlobalIndexer>
90 panzer::L2Projection::getTargetGlobalIndexer() const
91 {return targetGlobalIndexer_;}
92
93 Teuchos::RCP<Tpetra::CrsMatrix<double,panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType>>
94 panzer::L2Projection::buildMassMatrix(bool use_lumping,
95 const std::unordered_map<std::string,double>* elementBlockMultipliers)
96 {
97 PANZER_FUNC_TIME_MONITOR_DIFF("L2Projection::Build Mass Matrix",TopBuildMassMatrix);
98 TEUCHOS_ASSERT(setupCalled_);
99
100 if (elementBlockMultipliers != nullptr) {
101 TEUCHOS_ASSERT(elementBlockMultipliers->size() == elementBlockNames_.size());
102 }
103
104 // Allocate the owned matrix
105 std::vector<Teuchos::RCP<const panzer::GlobalIndexer>> indexers;
106 indexers.push_back(targetGlobalIndexer_);
107
108 panzer::BlockedTpetraLinearObjFactory<panzer::Traits,double,panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType> factory(comm_,indexers);
109
110 auto ownedMatrix = factory.getTpetraMatrix(0,0);
111 auto ghostedMatrix = factory.getGhostedTpetraMatrix(0,0);
112 ownedMatrix->resumeFill();
113 ownedMatrix->setAllToScalar(0.0);
114 ghostedMatrix->resumeFill();
115 ghostedMatrix->setAllToScalar(0.0);
116
117 const int fieldIndex = targetGlobalIndexer_->getFieldNum(targetBasisDescriptor_.getType());
118
119 const bool is_scalar = targetBasisDescriptor_.getType()=="HGrad" || targetBasisDescriptor_.getType()=="Const" || targetBasisDescriptor_.getType()=="HVol";
120
121 // Loop over element blocks and fill mass matrix
122 if(is_scalar){
123 auto M = ghostedMatrix->getLocalMatrixDevice();
124 for (const auto& block : elementBlockNames_) {
125
126 double ebMultiplier = 1.0;
127 if (elementBlockMultipliers != nullptr)
128 ebMultiplier = elementBlockMultipliers->find(block)->second;
129
130 // Get the cell local ids and set the BasisValues object (BV
131 // can come from a user defined set or from WorksetContainer).
132 const panzer::BasisValues2<double>* bv_ptr;
133 int num_cells_owned_ghosted_virtual = 0;
134 int num_cells_owned = 0;
135 Kokkos::View<const panzer::LocalOrdinal*> cell_local_ids;
136 if (useUserSuppliedBasisValues_) {
137 // Skip this block if there are no elements in this block on this mpi process
138 auto tmp = connManager_->getElementBlock(block); // no ghosting or virtual in this case
139 if (tmp.size() == 0)
140 continue;
141
142 Kokkos::View<panzer::LocalOrdinal*>::host_mirror_type cell_local_ids_host(tmp.data(),tmp.size());
143 Kokkos::View<panzer::LocalOrdinal*> cell_local_ids_nonconst("cell_local_ids",tmp.size());
144 Kokkos::deep_copy(cell_local_ids_nonconst,cell_local_ids_host);
145 cell_local_ids = cell_local_ids_nonconst;
146
147 bv_ptr = basisValues_[block].get();
148 num_cells_owned_ghosted_virtual = cell_local_ids.extent(0);
149 num_cells_owned = cell_local_ids.extent(0);
150 }
151 else {
152 // Based on descriptor, currently assumes there should only
153 // be one workset (partitioned path assumes a single
154 // workset).
155 panzer::WorksetDescriptor wd(block,panzer::WorksetSizeType::ALL_ELEMENTS,true,true);
156 const auto worksets = worksetContainer_->getWorksets(wd);
157 // Skip this block if there are no elements in this block on this mpi process
158 if (worksets->size() == 0)
159 continue;
160
161 TEUCHOS_ASSERT(worksets->size() == 1);
162
163 const auto& workset = (*worksets)[0];
164 bv_ptr = &(workset.getBasisValues(targetBasisDescriptor_,integrationDescriptor_));
165 num_cells_owned_ghosted_virtual = workset.numOwnedCells()+workset.numGhostCells()+workset.numVirtualCells();
166 num_cells_owned = workset.numOwnedCells();
167 cell_local_ids = workset.getLocalCellIDs();
168 }
169 const auto& basisValues = *bv_ptr;
170
171 const auto unweightedBasis = basisValues.getBasisValues(false).get_static_view();
172 const auto weightedBasis = basisValues.getBasisValues(true).get_static_view();
173
174 const std::vector<panzer::LocalOrdinal>& offsets = targetGlobalIndexer_->getGIDFieldOffsets(block,fieldIndex);
175 PHX::View<panzer::LocalOrdinal*> kOffsets("MassMatrix: Offsets",offsets.size());
176 auto kOffsets_h = Kokkos::create_mirror_view(kOffsets);
177
178 for(const auto& i : offsets)
179 kOffsets_h(i) = offsets[i];
180
181 Kokkos::deep_copy(kOffsets, kOffsets_h);
182
183 // Local Ids
184 PHX::View<panzer::LocalOrdinal**> localIds("MassMatrix: LocalIds", num_cells_owned_ghosted_virtual,
185 targetGlobalIndexer_->getElementBlockGIDCount(block));
186
187 // Remove the ghosted cell ids or the call to getElementLocalIds will spill array bounds
188 const auto cellLocalIdsNoGhost = Kokkos::subview(cell_local_ids,std::make_pair(0,num_cells_owned));
189
190 targetGlobalIndexer_->getElementLIDs(cellLocalIdsNoGhost,localIds);
191
192 const int numBasisPoints = static_cast<int>(weightedBasis.extent(1));
193 if ( use_lumping ) {
194 Kokkos::parallel_for(num_cells_owned,KOKKOS_LAMBDA (const int& cell) {
195 double total_mass = 0.0, trace = 0.0;
196
197 panzer::LocalOrdinal cLIDs[256];
198 const int numIds = static_cast<int>(localIds.extent(1));
199 for(int i=0;i<numIds;++i)
200 cLIDs[i] = localIds(cell,i);
201
202 double vals[256]={0.0};
203 const int numQP = static_cast<int>(unweightedBasis.extent(2));
204
205 for (int row=0; row < numBasisPoints; ++row) {
206 for (int col=0; col < numIds; ++col) {
207 for (int qp=0; qp < numQP; ++qp) {
208 auto tmp = unweightedBasis(cell,row,qp) * weightedBasis(cell,col,qp) * ebMultiplier;
209 total_mass += tmp;
210 if (col == row )
211 trace += tmp;
212 }
213 }
214 }
215
216 for (int row=0; row < numBasisPoints; ++row) {
217 for (int col=0; col < numBasisPoints; ++col)
218 vals[col] = 0;
219
220 int offset = kOffsets(row);
221 panzer::LocalOrdinal lid = localIds(cell,offset);
222 int col = row;
223 vals[col] = 0.0;
224 for (int qp=0; qp < numQP; ++qp)
225 vals[col] += unweightedBasis(cell,row,qp) * weightedBasis(cell,col,qp) * ebMultiplier * total_mass / trace;
226
227 M.sumIntoValues(lid,cLIDs,numIds,vals,true,true);
228 }
229 });
230
231 } else {
232 Kokkos::parallel_for(num_cells_owned,KOKKOS_LAMBDA (const int& cell) {
233 panzer::LocalOrdinal cLIDs[256];
234 const int numIds = static_cast<int>(localIds.extent(1));
235 for(int i=0;i<numIds;++i)
236 cLIDs[i] = localIds(cell,i);
237
238 double vals[256];
239 const int numQP = static_cast<int>(unweightedBasis.extent(2));
240
241 for (int row=0; row < numBasisPoints; ++row) {
242 int offset = kOffsets(row);
243 panzer::LocalOrdinal lid = localIds(cell,offset);
244
245 for (int col=0; col < numIds; ++col) {
246 vals[col] = 0.0;
247 for (int qp=0; qp < numQP; ++qp)
248 vals[col] += unweightedBasis(cell,row,qp) * weightedBasis(cell,col,qp) * ebMultiplier;
249 }
250 M.sumIntoValues(lid,cLIDs,numIds,vals,true,true);
251
252 }
253 });
254 }
255 }
256 } else {
257 auto M = ghostedMatrix->getLocalMatrixDevice();
258 for (const auto& block : elementBlockNames_) {
259
260 double ebMultiplier = 1.0;
261 if (elementBlockMultipliers != nullptr)
262 ebMultiplier = elementBlockMultipliers->find(block)->second;
263
264 // Get the cell local ids and set the BasisValues object (BV
265 // can come from a user defined set or from WorksetContainer).
266 const panzer::BasisValues2<double>* bv_ptr;
267 int num_cells_owned_ghosted_virtual = 0;
268 int num_cells_owned = 0;
269 Kokkos::View<const panzer::LocalOrdinal*> cell_local_ids;
270 if (useUserSuppliedBasisValues_) {
271 // Skip this block if there are no elements in this block on this mpi process
272 auto tmp = connManager_->getElementBlock(block); // no ghosting or virtual in this case
273 if (tmp.size() == 0)
274 continue;
275
276 Kokkos::View<panzer::LocalOrdinal*>::host_mirror_type cell_local_ids_host(tmp.data(),tmp.size());
277 Kokkos::View<panzer::LocalOrdinal*> cell_local_ids_nonconst("cell_local_ids",tmp.size());
278 Kokkos::deep_copy(cell_local_ids_nonconst,cell_local_ids_host);
279 cell_local_ids = cell_local_ids_nonconst;
280
281 bv_ptr = basisValues_[block].get();
282 num_cells_owned_ghosted_virtual = cell_local_ids.extent(0);
283 num_cells_owned = cell_local_ids.extent(0);
284 }
285 else {
286 // Based on descriptor, currently assumes there should only
287 // be one workset (partitioned path assumes a single
288 // workset).
289 panzer::WorksetDescriptor wd(block,panzer::WorksetSizeType::ALL_ELEMENTS,true,true);
290 const auto worksets = worksetContainer_->getWorksets(wd);
291 // Skip this block if there are no elements in this block on this mpi process
292 if (worksets->size() == 0)
293 continue;
294
295 TEUCHOS_ASSERT(worksets->size() == 1);
296
297 const auto& workset = (*worksets)[0];
298 bv_ptr = &(workset.getBasisValues(targetBasisDescriptor_,integrationDescriptor_));
299 num_cells_owned_ghosted_virtual = workset.numOwnedCells()+workset.numGhostCells()+workset.numVirtualCells();
300 num_cells_owned = workset.numOwnedCells();
301 cell_local_ids = workset.getLocalCellIDs();
302 }
303 const auto& basisValues = *bv_ptr;
304
305 const auto unweightedBasis = basisValues.getVectorBasisValues(false).get_static_view();
306 const auto weightedBasis = basisValues.getVectorBasisValues(true).get_static_view();
307
308 const std::vector<panzer::LocalOrdinal>& offsets = targetGlobalIndexer_->getGIDFieldOffsets(block,fieldIndex);
309 PHX::View<panzer::LocalOrdinal*> kOffsets("MassMatrix: Offsets",offsets.size());
310 auto kOffsets_h = Kokkos::create_mirror_view(kOffsets);
311
312 for(const auto& i : offsets)
313 kOffsets_h(i) = offsets[i];
314
315 Kokkos::deep_copy(kOffsets, kOffsets_h);
316
317 // Local Ids
318 PHX::View<panzer::LocalOrdinal**> localIds("MassMatrix: LocalIds", num_cells_owned_ghosted_virtual,
319 targetGlobalIndexer_->getElementBlockGIDCount(block));
320
321 // Remove the ghosted cell ids or the call to getElementLocalIds will spill array bounds
322 const PHX::View<const int*> cellLocalIdsNoGhost = Kokkos::subview(cell_local_ids,std::make_pair(0,num_cells_owned));
323
324 targetGlobalIndexer_->getElementLIDs(cellLocalIdsNoGhost,localIds);
325
326 const int numBasisPoints = static_cast<int>(weightedBasis.extent(1));
327 Kokkos::parallel_for(num_cells_owned,KOKKOS_LAMBDA (const int& cell) {
328
329 panzer::LocalOrdinal cLIDs[256];
330 const int numIds = static_cast<int>(localIds.extent(1));
331 for(int i=0;i<numIds;++i)
332 cLIDs[i] = localIds(cell,i);
333
334 double vals[256];
335 const int numQP = static_cast<int>(unweightedBasis.extent(2));
336
337 for (int qp=0; qp < numQP; ++qp) {
338 for (int row=0; row < numBasisPoints; ++row) {
339 int offset = kOffsets(row);
340 panzer::LocalOrdinal lid = localIds(cell,offset);
341
342 for (int col=0; col < numIds; ++col){
343 vals[col] = 0.0;
344 for(int dim=0; dim < static_cast<int>(weightedBasis.extent(3)); ++dim)
345 vals[col] += unweightedBasis(cell,row,qp,dim) * weightedBasis(cell,col,qp,dim) * ebMultiplier;
346 }
347
348 M.sumIntoValues(lid,cLIDs,numIds,vals,true,true);
349 }
350 }
351 });
352 }
353 }
354
355 {
356 PANZER_FUNC_TIME_MONITOR_DIFF("Exporting of mass matrix",ExportMM);
357 auto map = factory.getMap(0);
358 ghostedMatrix->fillComplete(map,map);
359 const auto exporter = factory.getGhostedExport(0);
360 ownedMatrix->doExport(*ghostedMatrix, *exporter, Tpetra::ADD);
361 ownedMatrix->fillComplete();
362 }
363 return ownedMatrix;
364 }
365
366 Teuchos::RCP<Tpetra::MultiVector<double,panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType>>
367 panzer::L2Projection::buildInverseLumpedMassMatrix()
368 {
369 PANZER_FUNC_TIME_MONITOR_DIFF("L2Projection<panzer::LocalOrdinal,panzer::GlobalOrdinal>::buildInverseLumpedMassMatrix",buildInvLMM);
370 using Teuchos::rcp;
371 const auto massMatrix = this->buildMassMatrix(true);
372 const auto lumpedMassMatrix = rcp(new Tpetra::MultiVector<double,panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType>(massMatrix->getDomainMap(),1,true));
373 const auto tmp = rcp(new Tpetra::MultiVector<double,panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType>(massMatrix->getRangeMap(),1,false));
374 tmp->putScalar(1.0);
375 {
376 PANZER_FUNC_TIME_MONITOR_DIFF("Apply",Apply);
377 massMatrix->apply(*tmp,*lumpedMassMatrix);
378 }
379 {
380 PANZER_FUNC_TIME_MONITOR_DIFF("Reciprocal",Reciprocal);
381 lumpedMassMatrix->reciprocal(*lumpedMassMatrix);
382 }
383 return lumpedMassMatrix;
384 }
385
386 Teuchos::RCP<Tpetra::CrsMatrix<double,panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType>>
387 panzer::L2Projection::buildRHSMatrix(const panzer::GlobalIndexer& sourceGlobalIndexer,
388 const Teuchos::RCP<const Tpetra::Map<panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType>>& inputOwnedSourceMap,
389 const std::string& sourceFieldName,
390 const panzer::BasisDescriptor& sourceBasisDescriptor,
391 const int directionIndex)
392 {
393 // *******************
394 // Create Retangular matrix (both ghosted and owned).
395 // *******************
396 using Teuchos::RCP;
397 using Teuchos::rcp;
398 using MapType = Tpetra::Map<panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType>;
399 using GraphType = Tpetra::CrsGraph<panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType>;
400 using ExportType = Tpetra::Export<panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType>;
401 using MatrixType = Tpetra::CrsMatrix<double,panzer::LocalOrdinal,panzer::GlobalOrdinal,panzer::TpetraNodeType>;
402
403 // *******************
404 // Build ghosted graph
405 // *******************
406
407 RCP<MapType> ghostedTargetMap;
408 {
409 std::vector<panzer::GlobalOrdinal> indices;
410 targetGlobalIndexer_->getOwnedAndGhostedIndices(indices);
411 ghostedTargetMap = rcp(new MapType(Teuchos::OrdinalTraits<panzer::GlobalOrdinal>::invalid(),indices,0,comm_));
412 }
413
414 RCP<MapType> ghostedSourceMap;
415 {
416 std::vector<panzer::GlobalOrdinal> indices;
417 sourceGlobalIndexer.getOwnedAndGhostedIndices(indices);
418 ghostedSourceMap = rcp(new MapType(Teuchos::OrdinalTraits<panzer::GlobalOrdinal>::invalid(),indices,0,comm_));
419 }
420
421
422 // Now insert the non-zero pattern per row
423 // count number of entries per row; required by CrsGraph constructor
424 std::vector<size_t> nEntriesPerRow(ghostedTargetMap->getLocalNumElements(),0);
425 std::vector<std::string> elementBlockIds;
426 targetGlobalIndexer_->getElementBlockIds(elementBlockIds);
427 std::vector<std::string>::const_iterator blockItr;
428 for (blockItr=elementBlockIds.begin();blockItr!=elementBlockIds.end();++blockItr) {
429 std::string blockId = *blockItr;
430 const std::vector<panzer::LocalOrdinal> & elements = targetGlobalIndexer_->getElementBlock(blockId);
431
432 std::vector<panzer::GlobalOrdinal> row_gids;
433 std::vector<panzer::GlobalOrdinal> col_gids;
434
435 for(std::size_t elmt=0;elmt<elements.size();elmt++) {
436 targetGlobalIndexer_->getElementGIDs(elements[elmt],row_gids);
437 sourceGlobalIndexer.getElementGIDs(elements[elmt],col_gids);
438 for(std::size_t row=0;row<row_gids.size();row++) {
439 panzer::LocalOrdinal lid =
440 ghostedTargetMap->getLocalElement(row_gids[row]);
441 nEntriesPerRow[lid] += col_gids.size();
442 }
443 }
444 }
445
446 Teuchos::ArrayView<const size_t> nEntriesPerRowView(nEntriesPerRow);
447 RCP<GraphType> ghostedGraph = rcp(new GraphType(ghostedTargetMap,ghostedSourceMap,nEntriesPerRowView));
448
449 for (blockItr=elementBlockIds.begin();blockItr!=elementBlockIds.end();++blockItr) {
450 std::string blockId = *blockItr;
451 const std::vector<panzer::LocalOrdinal> & elements = targetGlobalIndexer_->getElementBlock(blockId);
452
453 std::vector<panzer::GlobalOrdinal> row_gids;
454 std::vector<panzer::GlobalOrdinal> col_gids;
455
456 for(std::size_t elmt=0;elmt<elements.size();elmt++) {
457 targetGlobalIndexer_->getElementGIDs(elements[elmt],row_gids);
458 sourceGlobalIndexer.getElementGIDs(elements[elmt],col_gids);
459 for(std::size_t row=0;row<row_gids.size();row++)
460 ghostedGraph->insertGlobalIndices(row_gids[row],col_gids);
461 }
462 }
463
464 RCP<MapType> ownedTargetMap;
465 {
466 std::vector<panzer::GlobalOrdinal> indices;
467 targetGlobalIndexer_->getOwnedIndices(indices);
468 ownedTargetMap = rcp(new MapType(Teuchos::OrdinalTraits<panzer::GlobalOrdinal>::invalid(),indices,0,comm_));
469 }
470
471 RCP<const MapType> ownedSourceMap = inputOwnedSourceMap;
472 if (is_null(ownedSourceMap)) {
473 std::vector<panzer::GlobalOrdinal> indices;
474 sourceGlobalIndexer.getOwnedIndices(indices);
475 ownedSourceMap = rcp(new MapType(Teuchos::OrdinalTraits<panzer::GlobalOrdinal>::invalid(),indices,0,comm_));
476 }
477
478 // Fill complete with owned range and domain map
479 ghostedGraph->fillComplete(ownedSourceMap,ownedTargetMap);
480
481 // *****************
482 // Build owned graph
483 // *****************
484
485 RCP<GraphType> ownedGraph = rcp(new GraphType(ownedTargetMap,0));
486 RCP<const ExportType> exporter = rcp(new ExportType(ghostedTargetMap,ownedTargetMap));
487 ownedGraph->doExport(*ghostedGraph, *exporter, Tpetra::INSERT);
488 ownedGraph->fillComplete(ownedSourceMap,ownedTargetMap);
489
490 RCP<MatrixType> ghostedMatrix = rcp(new MatrixType(ghostedGraph));
491 RCP<MatrixType> ownedMatrix = rcp(new MatrixType(ownedGraph));
492 // ghostedMatrix.fillComplete();
493 // ghostedMatrix.resumeFill();
494
495 ghostedMatrix->setAllToScalar(0.0);
496 ownedMatrix->setAllToScalar(0.0);
497
498 // *******************
499 // Fill ghosted matrix
500 // *******************
501 for (const auto& block : elementBlockNames_) {
502
503 panzer::WorksetDescriptor wd(block,panzer::WorksetSizeType::ALL_ELEMENTS,true,true);
504 const auto& worksets = worksetContainer_->getWorksets(wd);
505 for (const auto& workset : *worksets) {
506
507 // Get target basis values: current implementation assumes target basis is HGrad
508 const auto& targetBasisValues = workset.getBasisValues(targetBasisDescriptor_,integrationDescriptor_);
509 const auto& targetWeightedBasis = targetBasisValues.getBasisValues(true).get_static_view();
510
511 // Sources can be any basis
512 const auto& sourceBasisValues = workset.getBasisValues(sourceBasisDescriptor,integrationDescriptor_);
513 PHX::View<const double***> sourceUnweightedScalarBasis;
514 PHX::View<const double****> sourceUnweightedVectorBasis;
515 bool useRankThreeBasis = false; // default to gradient or vector basis
516 if ( (sourceBasisDescriptor.getType() == "HGrad") || (sourceBasisDescriptor.getType() == "Const") || (sourceBasisDescriptor.getType() == "HVol") ) {
517 if (directionIndex == -1) { // Project dof value
518 sourceUnweightedScalarBasis = sourceBasisValues.getBasisValues(false).get_static_view();
519 useRankThreeBasis = true;
520 }
521 else { // Project dof gradient of scalar basis
522 sourceUnweightedVectorBasis = sourceBasisValues.getGradBasisValues(false).get_static_view();
523 }
524 }
525 else { // Project vector value
526 sourceUnweightedVectorBasis = sourceBasisValues.getVectorBasisValues(false).get_static_view();
527 }
528
529 // Get the element local ids
530 PHX::View<panzer::LocalOrdinal**> targetLocalIds("buildRHSMatrix: targetLocalIds", workset.numOwnedCells(),
531 targetGlobalIndexer_->getElementBlockGIDCount(block));
532 PHX::View<panzer::LocalOrdinal**> sourceLocalIds("buildRHSMatrix: sourceLocalIds", workset.numOwnedCells(),
533 sourceGlobalIndexer.getElementBlockGIDCount(block));
534 {
535 // Remove the ghosted cell ids or the call to getElementLocalIds will spill array bounds
536 const auto cellLocalIdsNoGhost = Kokkos::subview(workset.cell_local_ids_k,std::make_pair(0,workset.numOwnedCells()));
537 targetGlobalIndexer_->getElementLIDs(cellLocalIdsNoGhost,targetLocalIds);
538 sourceGlobalIndexer.getElementLIDs(cellLocalIdsNoGhost,sourceLocalIds);
539 }
540
541 // Get the offsets
542 PHX::View<panzer::LocalOrdinal*> targetFieldOffsets;
543 {
544 const auto fieldIndex = targetGlobalIndexer_->getFieldNum(targetBasisDescriptor_.getType());
545 const std::vector<panzer::LocalOrdinal>& offsets = targetGlobalIndexer_->getGIDFieldOffsets(block,fieldIndex);
546 targetFieldOffsets = PHX::View<panzer::LocalOrdinal*>("L2Projection:buildRHS:targetFieldOffsets",offsets.size());
547 const auto hostOffsets = Kokkos::create_mirror_view(targetFieldOffsets);
548 for(size_t i=0; i < offsets.size(); ++i)
549 hostOffsets(i) = offsets[i];
550 Kokkos::deep_copy(targetFieldOffsets,hostOffsets);
551 }
552
553 PHX::View<panzer::LocalOrdinal*> sourceFieldOffsets;
554 {
555 const auto fieldIndex = sourceGlobalIndexer.getFieldNum(sourceFieldName);
556 const std::vector<panzer::LocalOrdinal>& offsets = sourceGlobalIndexer.getGIDFieldOffsets(block,fieldIndex);
557 TEUCHOS_TEST_FOR_EXCEPTION(offsets.size() == 0, std::runtime_error,
558 "ERROR: panzer::L2Projection::buildRHSMatrix() - The source field, \""
559 << sourceFieldName << "\", does not exist in element block \""
560 << block << "\"!");
561 sourceFieldOffsets = PHX::View<panzer::LocalOrdinal*>("L2Projection:buildRHS:sourceFieldOffsets",offsets.size());
562 const auto hostOffsets = Kokkos::create_mirror_view(sourceFieldOffsets);
563 for(size_t i=0; i <offsets.size(); ++i)
564 hostOffsets(i) = offsets[i];
565 Kokkos::deep_copy(sourceFieldOffsets,hostOffsets);
566 }
567
568 const auto localMatrix = ghostedMatrix->getLocalMatrixDevice();
569 const int numRows = static_cast<int>(targetWeightedBasis.extent(1));
570 int tmpNumCols = -1;
571 int tmpNumQP = -1;
572 if (useRankThreeBasis) {
573 tmpNumCols = static_cast<int>(sourceUnweightedScalarBasis.extent(1));
574 tmpNumQP = static_cast<int>(sourceUnweightedScalarBasis.extent(2));
575 }
576 else {
577 tmpNumCols = static_cast<int>(sourceUnweightedVectorBasis.extent(1));
578 tmpNumQP = static_cast<int>(sourceUnweightedVectorBasis.extent(2));
579 }
580 const int numCols = tmpNumCols;
581 const int numQP = tmpNumQP;
582
583 Kokkos::parallel_for(Kokkos::RangePolicy<PHX::Device>(0,workset.numOwnedCells()),KOKKOS_LAMBDA (const int& cell) {
584 panzer::LocalOrdinal cLIDs[256];
585 double vals[256];
586 for (int row = 0; row < numRows; ++row) {
587 const int rowOffset = targetFieldOffsets(row);
588 const int rowLID = targetLocalIds(cell,rowOffset);
589 for (int col = 0; col < numCols; ++col)
590 vals[col] = 0.0;
591
592 for (int col = 0; col < numCols; ++col) {
593 for (int qp = 0; qp < numQP; ++qp) {
594 const int colOffset = sourceFieldOffsets(col);
595 const int colLID = sourceLocalIds(cell,colOffset);
596 cLIDs[col] = colLID;
597 if (useRankThreeBasis)
598 vals[col] += sourceUnweightedScalarBasis(cell,col,qp) * targetWeightedBasis(cell,row,qp);
599 else
600 vals[col] += sourceUnweightedVectorBasis(cell,col,qp,directionIndex) * targetWeightedBasis(cell,row,qp);
601 }
602 }
603 localMatrix.sumIntoValues(rowLID,cLIDs,numCols,vals,true,true);
604 }
605 });
606 }
607
608 }
609 ghostedMatrix->fillComplete(ownedSourceMap,ownedTargetMap);
610 ownedMatrix->resumeFill();
611 ownedMatrix->doExport(*ghostedMatrix,*exporter,Tpetra::ADD);
612 ownedMatrix->fillComplete(ownedSourceMap,ownedTargetMap);
613
614 return ownedMatrix;
615 }
616
617}
618
619#endif
Panzer_BasisDescriptor.hpp
Panzer_BlockedDOFManagerFactory.hpp
Panzer_BlockedTpetraLinearObjFactory.hpp
offsets
PHX::View< const int * > offsets
Definition Panzer_DOFCurl_impl.hpp:148
Panzer_DOFManagerFactory.hpp
Panzer_ElementBlockIdToPhysicsIdMap.hpp
Panzer_IntrepidFieldPattern.hpp
Panzer_L2Projection.hpp
vals
Kokkos::View< typename Sacado::ScalarType< ScalarT >::type **, Kokkos::LayoutRight, PHX::Device > vals
Definition Panzer_ScatterResidual_Tpetra_impl.hpp:340
Panzer_TpetraLinearObjContainer.hpp
Panzer_WorksetContainer.hpp
Panzer_WorksetDescriptor.hpp
Panzer_Workset.hpp
panzer::BasisDescriptor::getType
const std::string & getType() const
Get type of basis.
Definition Panzer_BasisDescriptor.hpp:46
panzer::BasisValues2::getBasisValues
ConstArray_CellBasisIP getBasisValues(const bool weighted, const bool cache=true, const bool force=false) const
Get the basis values evaluated at mesh points.
Definition Panzer_BasisValues2_impl.hpp:1035
panzer::BasisValues2::getVectorBasisValues
ConstArray_CellBasisIPDim getVectorBasisValues(const bool weighted, const bool cache=true, const bool force=false) const
Get the vector basis values evaluated at mesh points.
Definition Panzer_BasisValues2_impl.hpp:1202
panzer::BlockedTpetraLinearObjFactory::getMap
virtual Teuchos::RCP< const MapType > getMap(int i) const
get the map from the matrix
Definition Panzer_BlockedTpetraLinearObjFactory_impl.hpp:829
panzer::BlockedTpetraLinearObjFactory::getTpetraMatrix
Teuchos::RCP< CrsMatrixType > getTpetraMatrix(int i, int j) const
Definition Panzer_BlockedTpetraLinearObjFactory_impl.hpp:1047
panzer::BlockedTpetraLinearObjFactory::getGhostedTpetraMatrix
Teuchos::RCP< CrsMatrixType > getGhostedTpetraMatrix(int i, int j) const
Definition Panzer_BlockedTpetraLinearObjFactory_impl.hpp:1062
panzer::BlockedTpetraLinearObjFactory::getGhostedExport
virtual Teuchos::RCP< const ExportType > getGhostedExport(int j) const
get exporter for converting an overalapped object to a "normal" object
Definition Panzer_BlockedTpetraLinearObjFactory_impl.hpp:903
panzer::GlobalIndexer::getGIDFieldOffsets
virtual const std::vector< int > & getGIDFieldOffsets(const std::string &blockId, int fieldNum) const =0
Use the field pattern so that you can find a particular field in the GIDs array.
panzer::GlobalIndexer::getElementGIDs
virtual void getElementGIDs(panzer::LocalOrdinal localElmtId, std::vector< panzer::GlobalOrdinal > &gids, const std::string &blockIdHint="") const =0
Get the global IDs for a particular element. This function overwrites the gids variable.
panzer::GlobalIndexer::getOwnedAndGhostedIndices
virtual void getOwnedAndGhostedIndices(std::vector< panzer::GlobalOrdinal > &indices) const =0
Get the set of owned and ghosted indices for this processor.
panzer::GlobalIndexer::getElementBlockGIDCount
virtual int getElementBlockGIDCount(const std::size_t &blockIndex) const =0
How any GIDs are associate with each element in a particular element block.
panzer::GlobalIndexer::getOwnedIndices
virtual void getOwnedIndices(std::vector< panzer::GlobalOrdinal > &indices) const =0
Get the set of indices owned by this processor.
panzer::GlobalIndexer::getElementLIDs
const Kokkos::View< const panzer::LocalOrdinal *, Kokkos::LayoutRight, PHX::Device > getElementLIDs(panzer::LocalOrdinal localElmtId) const
Definition Panzer_GlobalIndexer.hpp:221
panzer::GlobalIndexer::getFieldNum
virtual int getFieldNum(const std::string &str) const =0
Get the number used for access to this field.
panzer::L2Projection::useBasisValues
void useBasisValues(const std::map< std::string, Teuchos::RCP< panzer::BasisValues2< double > > > &map_eblock_to_bv)
Override using the panzer::WorksetContainer and instead use the registered BasisValues object.
Definition Panzer_L2Projection.cpp:77
panzer::L2Projection::buildMassMatrix
Teuchos::RCP< Tpetra::CrsMatrix< double, panzer::LocalOrdinal, panzer::GlobalOrdinal, Tpetra::KokkosCompat::KokkosDeviceWrapperNode< PHX::Device > > > buildMassMatrix(bool use_lumping=false, const std::unordered_map< std::string, double > *elementBlockMultipliers=nullptr)
Allocates, fills and returns a mass matrix for L2 projection onto a target basis.
Definition Panzer_L2Projection.cpp:94
panzer::L2Projection::getTargetGlobalIndexer
Teuchos::RCP< panzer::GlobalIndexer > getTargetGlobalIndexer() const
Returns the target global indexer. Will be null if setup() has not been called.
Definition Panzer_L2Projection.cpp:90
panzer::L2Projection::buildInverseLumpedMassMatrix
Teuchos::RCP< Tpetra::MultiVector< double, panzer::LocalOrdinal, panzer::GlobalOrdinal, Tpetra::KokkosCompat::KokkosDeviceWrapperNode< PHX::Device > > > buildInverseLumpedMassMatrix()
Allocates, fills and returns a Tpetra::MultiVector containing the inverse lumped mass matrix values a...
Definition Panzer_L2Projection.cpp:367
panzer::L2Projection::buildRHSMatrix
Teuchos::RCP< Tpetra::CrsMatrix< double, panzer::LocalOrdinal, panzer::GlobalOrdinal, Tpetra::KokkosCompat::KokkosDeviceWrapperNode< PHX::Device > > > buildRHSMatrix(const panzer::GlobalIndexer &sourceDOFManager, const Teuchos::RCP< const Tpetra::Map< panzer::LocalOrdinal, panzer::GlobalOrdinal, Tpetra::KokkosCompat::KokkosDeviceWrapperNode< PHX::Device > > > &ownedSourceMap, const std::string &sourceFieldName, const panzer::BasisDescriptor &sourceBasisDescriptor, const int vectorOrGradientDirectionIndex=-1)
Allocates, fills and returns a rectangular matrix for L2 projection of a scalar field,...
Definition Panzer_L2Projection.cpp:387
panzer::L2Projection::setup
void setup(const panzer::BasisDescriptor &targetBasis, const panzer::IntegrationDescriptor &integrationDescriptor, const Teuchos::RCP< const Teuchos::MpiComm< int > > &comm, const Teuchos::RCP< const panzer::ConnManager > &connManager, const std::vector< std::string > &elementBlockNames, const Teuchos::RCP< panzer::WorksetContainer > worksetContainer=Teuchos::null)
Setup base objects for L2 Projections - requires target scalar basis and creates worksets if not supp...
Definition Panzer_L2Projection.cpp:33
panzer::ALL_ELEMENTS
@ ALL_ELEMENTS
Workset size is set to the total number of local elements in the MPI process.
Definition Panzer_WorksetDescriptor.hpp:26
Generated by doxygen 1.9.8