4 #ifndef PANZER_L2_PROJECTION_IMPL_HPP 5 #define PANZER_L2_PROJECTION_IMPL_HPP 7 #include "Teuchos_DefaultMpiComm.hpp" 8 #include "Tpetra_CrsGraph.hpp" 9 #include "Tpetra_MultiVector.hpp" 10 #include "Tpetra_CrsMatrix.hpp" 13 #include "Panzer_TpetraLinearObjFactory.hpp" 16 #include "Panzer_DOFManagerFactory.hpp" 17 #include "Panzer_BlockedDOFManagerFactory.hpp" 25 template<
typename LO,
typename GO>
30 const std::vector<std::string>& elementBlockNames,
31 const Teuchos::RCP<panzer::WorksetContainer> worksetContainer)
33 targetBasisDescriptor_ = targetBasis;
34 integrationDescriptor_ = integrationDescriptor;
36 connManager_ = connManager;
37 elementBlockNames_ = elementBlockNames;
38 worksetContainer_ = worksetContainer;
42 targetGlobalIndexer_ =
46 for (
const auto& eBlock : elementBlockNames_) {
47 std::vector<shards::CellTopology> topologies;
48 connManager_->getElementBlockTopologies(topologies);
49 std::vector<std::string> ebNames;
50 connManager_->getElementBlockIds(ebNames);
51 const auto search = std::find(ebNames.cbegin(),ebNames.cend(),eBlock);
52 TEUCHOS_ASSERT(search != ebNames.cend());
53 const int index = std::distance(ebNames.cbegin(),search);
54 const auto& cellTopology = topologies[index];
56 auto intrepidBasis = panzer::createIntrepid2Basis<PHX::Device,double,double>(targetBasisDescriptor_.getType(),
57 targetBasisDescriptor_.getOrder(),
61 targetGlobalIndexer_->addField(eBlock,targetBasisDescriptor_.getType(),fieldPattern);
64 targetGlobalIndexer_->buildGlobalUnknowns();
69 template<
typename LO,
typename GO>
70 Teuchos::RCP<panzer::UniqueGlobalIndexer<LO,GO>>
72 {
return targetGlobalIndexer_;}
74 template<
typename LO,
typename GO>
75 Teuchos::RCP<Tpetra::CrsMatrix<double,LO,GO,Kokkos::Compat::KokkosDeviceWrapperNode<PHX::Device>>>
78 TEUCHOS_ASSERT(setupCalled_);
81 std::vector<Teuchos::RCP<const panzer::UniqueGlobalIndexer<LO,GO>>> indexers;
82 indexers.push_back(targetGlobalIndexer_);
88 ownedMatrix->resumeFill();
89 ownedMatrix->setAllToScalar(0.0);
90 ghostedMatrix->resumeFill();
91 ghostedMatrix->setAllToScalar(0.0);
94 auto M = ghostedMatrix->getLocalMatrix();
95 const int fieldIndex = targetGlobalIndexer_->getFieldNum(targetBasisDescriptor_.getType());
97 const bool is_scalar = targetBasisDescriptor_.getType()==
"HGrad" || targetBasisDescriptor_.getType()==
"Const";
101 for (
const auto& block : elementBlockNames_) {
105 const auto& worksets = worksetContainer_->getWorksets(wd);
107 for (
const auto& workset : *worksets) {
109 const auto& basisValues = workset.getBasisValues(targetBasisDescriptor_,integrationDescriptor_);
111 const auto& unweightedBasis = basisValues.basis_scalar;
112 const auto& weightedBasis = basisValues.weighted_basis_scalar;
115 const std::vector<LO>&
offsets = targetGlobalIndexer_->getGIDFieldOffsets(block,fieldIndex);
116 PHX::View<LO*> kOffsets(
"MassMatrix: Offsets",
offsets.size());
120 PHX::Device::fence();
123 Kokkos::View<LO**,PHX::Device> localIds(
"MassMatrix: LocalIds", workset.numOwnedCells()+workset.numGhostCells()+workset.numVirtualCells(),
124 targetGlobalIndexer_->getElementBlockGIDCount(block));
127 const auto cellLocalIdsNoGhost = Kokkos::subview(workset.cell_local_ids_k,std::make_pair(0,workset.numOwnedCells()));
129 targetGlobalIndexer_->getElementLIDs(cellLocalIdsNoGhost,localIds);
131 const int numBasisPoints =
static_cast<int>(weightedBasis.extent(1));
132 Kokkos::parallel_for(workset.numOwnedCells(),KOKKOS_LAMBDA (
const int& cell) {
135 const int numIds =
static_cast<int>(localIds.extent(1));
136 for(
int i=0;i<numIds;++i)
137 cLIDs[i] = localIds(cell,i);
140 const int numQP =
static_cast<int>(unweightedBasis.extent(2));
142 for (
int qp=0; qp < numQP; ++qp) {
143 for (
int row=0; row < numBasisPoints; ++row) {
144 int offset = kOffsets(row);
145 LO lid = localIds(cell,offset);
147 for (
int col=0; col < numIds; ++col)
148 vals[col] = unweightedBasis(cell,row,qp) * weightedBasis(cell,col,qp);
150 M.sumIntoValues(lid,cLIDs,numIds,vals,
true,
true);
158 for (
const auto& block : elementBlockNames_) {
162 const auto& worksets = worksetContainer_->getWorksets(wd);
164 for (
const auto& workset : *worksets) {
166 const auto& basisValues = workset.getBasisValues(targetBasisDescriptor_,integrationDescriptor_);
168 const auto& unweightedBasis = basisValues.basis_vector;
169 const auto& weightedBasis = basisValues.weighted_basis_vector;
172 const std::vector<LO>&
offsets = targetGlobalIndexer_->getGIDFieldOffsets(block,fieldIndex);
173 PHX::View<LO*> kOffsets(
"MassMatrix: Offsets",
offsets.size());
177 PHX::Device::fence();
180 Kokkos::View<LO**,PHX::Device> localIds(
"MassMatrix: LocalIds", workset.numOwnedCells()+workset.numGhostCells()+workset.numVirtualCells(),
181 targetGlobalIndexer_->getElementBlockGIDCount(block));
184 const auto cellLocalIdsNoGhost = Kokkos::subview(workset.cell_local_ids_k,std::make_pair(0,workset.numOwnedCells()));
186 targetGlobalIndexer_->getElementLIDs(cellLocalIdsNoGhost,localIds);
188 const int numBasisPoints =
static_cast<int>(weightedBasis.extent(1));
189 Kokkos::parallel_for(workset.numOwnedCells(),KOKKOS_LAMBDA (
const int& cell) {
192 const int numIds =
static_cast<int>(localIds.extent(1));
193 for(
int i=0;i<numIds;++i)
194 cLIDs[i] = localIds(cell,i);
197 const int numQP =
static_cast<int>(unweightedBasis.extent(2));
199 for (
int qp=0; qp < numQP; ++qp) {
200 for (
int row=0; row < numBasisPoints; ++row) {
201 int offset = kOffsets(row);
202 LO lid = localIds(cell,offset);
204 for (
int col=0; col < numIds; ++col){
206 for(
int dim=0; dim < weightedBasis.extent(3); dim++)
207 vals[col] += unweightedBasis(cell,row,qp,dim) * weightedBasis(cell,col,qp,dim);
210 M.sumIntoValues(lid,cLIDs,numIds,vals,
true,
true);
218 PHX::exec_space::fence();
220 ghostedMatrix->fillComplete();
222 ownedMatrix->doExport(*ghostedMatrix, *exporter, Tpetra::ADD);
223 ownedMatrix->fillComplete();
228 template<
typename LO,
typename GO>
229 Teuchos::RCP<Tpetra::MultiVector<double,LO,GO,Kokkos::Compat::KokkosDeviceWrapperNode<PHX::Device>>>
233 const auto massMatrix = this->buildMassMatrix();
234 const auto lumpedMassMatrix = rcp(
new Tpetra::MultiVector<
double,LO,GO,Kokkos::Compat::KokkosDeviceWrapperNode<PHX::Device>>(massMatrix->getDomainMap(),1,
true));
235 const auto tmp = rcp(
new Tpetra::MultiVector<
double,LO,GO,Kokkos::Compat::KokkosDeviceWrapperNode<PHX::Device>>(massMatrix->getRangeMap(),1,
false));
237 massMatrix->apply(*tmp,*lumpedMassMatrix);
238 lumpedMassMatrix->reciprocal(*lumpedMassMatrix);
239 return lumpedMassMatrix;
242 template<
typename LO,
typename GO>
243 Teuchos::RCP<Tpetra::CrsMatrix<double,LO,GO,Kokkos::Compat::KokkosDeviceWrapperNode<PHX::Device>>>
245 const Teuchos::RCP<
const Tpetra::Map<LO,GO,Kokkos::Compat::KokkosDeviceWrapperNode<PHX::Device>>>& inputOwnedSourceMap,
246 const std::string& sourceFieldName,
248 const int directionIndex)
255 using MapType = Tpetra::Map<LO,GO,Kokkos::Compat::KokkosDeviceWrapperNode<PHX::Device>>;
256 using GraphType = Tpetra::CrsGraph<LO,GO,Kokkos::Compat::KokkosDeviceWrapperNode<PHX::Device>>;
257 using ExportType = Tpetra::Export<LO,GO,Kokkos::Compat::KokkosDeviceWrapperNode<PHX::Device>>;
258 using MatrixType = Tpetra::CrsMatrix<double,LO,GO,Kokkos::Compat::KokkosDeviceWrapperNode<PHX::Device>>;
264 RCP<MapType> ghostedTargetMap;
266 std::vector<GO> indices;
267 targetGlobalIndexer_->getOwnedAndGhostedIndices(indices);
268 ghostedTargetMap = rcp(
new MapType(Teuchos::OrdinalTraits<GO>::invalid(),indices,0,comm_));
271 RCP<MapType> ghostedSourceMap;
273 std::vector<GO> indices;
275 ghostedSourceMap = rcp(
new MapType(Teuchos::OrdinalTraits<GO>::invalid(),indices,0,comm_));
278 RCP<GraphType> ghostedGraph = rcp(
new GraphType(ghostedTargetMap,ghostedSourceMap,0));
281 std::vector<std::string> elementBlockIds;
282 targetGlobalIndexer_->getElementBlockIds(elementBlockIds);
283 std::vector<std::string>::const_iterator blockItr;
284 for (blockItr=elementBlockIds.begin();blockItr!=elementBlockIds.end();++blockItr) {
285 std::string blockId = *blockItr;
286 const std::vector<LO> & elements = targetGlobalIndexer_->getElementBlock(blockId);
288 std::vector<GO> row_gids;
289 std::vector<GO> col_gids;
291 for(std::size_t elmt=0;elmt<elements.size();elmt++) {
292 targetGlobalIndexer_->getElementGIDs(elements[elmt],row_gids);
294 for(std::size_t row=0;row<row_gids.size();row++)
295 ghostedGraph->insertGlobalIndices(row_gids[row],col_gids);
300 ghostedGraph->fillComplete(ghostedSourceMap,ghostedTargetMap);
305 RCP<MapType> ownedTargetMap;
307 std::vector<GO> indices;
308 targetGlobalIndexer_->getOwnedIndices(indices);
309 ownedTargetMap = rcp(
new MapType(Teuchos::OrdinalTraits<GO>::invalid(),indices,0,comm_));
312 RCP<const MapType> ownedSourceMap = inputOwnedSourceMap;
313 if (is_null(ownedSourceMap)) {
314 std::vector<GO> indices;
316 ownedSourceMap = rcp(
new MapType(Teuchos::OrdinalTraits<GO>::invalid(),indices,0,comm_));
319 RCP<GraphType> ownedGraph = rcp(
new GraphType(ownedTargetMap,0));
320 RCP<const ExportType> exporter = rcp(
new ExportType(ghostedTargetMap,ownedTargetMap));
321 ownedGraph->doExport(*ghostedGraph, *exporter, Tpetra::INSERT);
322 ownedGraph->fillComplete(ownedSourceMap,ownedTargetMap);
324 RCP<MatrixType> ghostedMatrix = rcp(
new MatrixType(ghostedGraph));
325 RCP<MatrixType> ownedMatrix = rcp(
new MatrixType(ownedGraph));
329 ghostedMatrix->setAllToScalar(0.0);
330 ownedMatrix->setAllToScalar(0.0);
331 PHX::Device::fence();
336 for (
const auto& block : elementBlockNames_) {
339 const auto& worksets = worksetContainer_->getWorksets(wd);
340 for (
const auto& workset : *worksets) {
343 const auto& targetBasisValues = workset.getBasisValues(targetBasisDescriptor_,integrationDescriptor_);
344 const auto& targetWeightedBasis = targetBasisValues.weighted_basis_scalar.get_static_view();
347 const auto& sourceBasisValues = workset.getBasisValues(sourceBasisDescriptor,integrationDescriptor_);
348 Kokkos::View<const double***,PHX::Device> sourceUnweightedScalarBasis;
349 Kokkos::View<const double****,PHX::Device> sourceUnweightedVectorBasis;
350 bool useRankThreeBasis =
false;
351 if ( (sourceBasisDescriptor.
getType() ==
"HGrad") || (sourceBasisDescriptor.
getType() ==
"Const") ) {
352 if (directionIndex == -1) {
353 sourceUnweightedScalarBasis = sourceBasisValues.basis_scalar.get_static_view();
354 useRankThreeBasis =
true;
357 sourceUnweightedVectorBasis = sourceBasisValues.grad_basis.get_static_view();
361 sourceUnweightedVectorBasis = sourceBasisValues.basis_vector.get_static_view();
365 Kokkos::View<LO**,PHX::Device> targetLocalIds(
"buildRHSMatrix: targetLocalIds", workset.numOwnedCells(),
366 targetGlobalIndexer_->getElementBlockGIDCount(block));
367 Kokkos::View<LO**,PHX::Device> sourceLocalIds(
"buildRHSMatrix: sourceLocalIds", workset.numOwnedCells(),
371 const auto cellLocalIdsNoGhost = Kokkos::subview(workset.cell_local_ids_k,std::make_pair(0,workset.numOwnedCells()));
372 targetGlobalIndexer_->getElementLIDs(cellLocalIdsNoGhost,targetLocalIds);
373 sourceGlobalIndexer.
getElementLIDs(cellLocalIdsNoGhost,sourceLocalIds);
377 Kokkos::View<LO*,PHX::Device> targetFieldOffsets;
379 const auto fieldIndex = targetGlobalIndexer_->getFieldNum(targetBasisDescriptor_.getType());
380 const std::vector<LO>&
offsets = targetGlobalIndexer_->getGIDFieldOffsets(block,fieldIndex);
381 targetFieldOffsets = Kokkos::View<LO*,PHX::Device>(
"L2Projection:buildRHS:targetFieldOffsets",
offsets.size());
382 const auto hostOffsets = Kokkos::create_mirror_view(targetFieldOffsets);
383 for(
size_t i=0; i <
offsets.size(); ++i)
385 Kokkos::deep_copy(targetFieldOffsets,hostOffsets);
386 PHX::Device::fence();
389 Kokkos::View<LO*,PHX::Device> sourceFieldOffsets;
391 const auto fieldIndex = sourceGlobalIndexer.
getFieldNum(sourceFieldName);
393 sourceFieldOffsets = Kokkos::View<LO*,PHX::Device>(
"L2Projection:buildRHS:sourceFieldOffsets",
offsets.size());
394 const auto hostOffsets = Kokkos::create_mirror_view(sourceFieldOffsets);
395 for(
size_t i=0; i <
offsets.size(); ++i)
397 Kokkos::deep_copy(sourceFieldOffsets,hostOffsets);
398 PHX::Device::fence();
401 const auto localMatrix = ghostedMatrix->getLocalMatrix();
402 const int numRows =
static_cast<int>(targetWeightedBasis.extent(1));
405 if (useRankThreeBasis) {
406 tmpNumCols =
static_cast<int>(sourceUnweightedScalarBasis.extent(1));
407 tmpNumQP =
static_cast<int>(sourceUnweightedScalarBasis.extent(2));
410 tmpNumCols =
static_cast<int>(sourceUnweightedVectorBasis.extent(1));
411 tmpNumQP =
static_cast<int>(sourceUnweightedVectorBasis.extent(2));
413 const int numCols = tmpNumCols;
414 const int numQP = tmpNumQP;
416 Kokkos::parallel_for(Kokkos::RangePolicy<PHX::Device>(0,workset.numOwnedCells()),KOKKOS_LAMBDA (
const int& cell) {
419 for (
int row = 0; row < numRows; ++row) {
420 const int rowOffset = targetFieldOffsets(row);
421 const int rowLID = targetLocalIds(cell,rowOffset);
422 for (
int col = 0; col < numCols; ++col)
425 for (
int col = 0; col < numCols; ++col) {
426 for (
int qp = 0; qp < numQP; ++qp) {
427 const int colOffset = sourceFieldOffsets(col);
428 const int colLID = sourceLocalIds(cell,colOffset);
430 if (useRankThreeBasis)
431 vals[col] += sourceUnweightedScalarBasis(cell,col,qp) * targetWeightedBasis(cell,row,qp);
433 vals[col] += sourceUnweightedVectorBasis(cell,col,qp,directionIndex) * targetWeightedBasis(cell,row,qp);
436 localMatrix.sumIntoValues(rowLID,cLIDs,numCols,vals,
true,
true);
442 ghostedMatrix->fillComplete(ghostedSourceMap,ghostedTargetMap);
443 ownedMatrix->resumeFill();
444 ownedMatrix->doExport(*ghostedMatrix,*exporter,Tpetra::ADD);
445 ownedMatrix->fillComplete(ownedSourceMap,ownedTargetMap);
virtual void getElementGIDs(LocalOrdinalT localElmtId, std::vector< GlobalOrdinalT > &gids, const std::string &blockIdHint="") const =0
Get the global IDs for a particular element. This function overwrites the gids variable.
virtual void getOwnedAndGhostedIndices(std::vector< GlobalOrdinalT > &indices) const =0
Get the set of owned and ghosted indices for this processor.
virtual Teuchos::RCP< const ExportType > getGhostedExport(int j) const
get exporter for converting an overalapped object to a "normal" object
const Kokkos::View< const LocalOrdinalT *, Kokkos::LayoutRight, PHX::Device > getElementLIDs(LocalOrdinalT localElmtId) const
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.
virtual int getElementBlockGIDCount(const std::string &blockId) const =0
How many GIDs are associate with a particular element block.
Teuchos::RCP< panzer::UniqueGlobalIndexer< LO, GO > > getTargetGlobalIndexer() const
Returns the target global indexer. Will be null if setup() has not been called.
virtual int getFieldNum(const std::string &str) const =0
Get the number used for access to this field.
Teuchos::RCP< Tpetra::MultiVector< double, LO, GO, Kokkos::Compat::KokkosDeviceWrapperNode< PHX::Device > > > buildInverseLumpedMassMatrix()
Allocates, fills and returns a Tpetra::MultiVector containing the inverse lumped mass matrix values...
virtual void getOwnedIndices(std::vector< GlobalOrdinalT > &indices) const =0
Get the set of indices owned by this processor.
void setup(const panzer::BasisDescriptor &targetBasis, const panzer::IntegrationDescriptor &integrationDescriptor, const Teuchos::RCP< const Teuchos::MpiComm< int >> &comm, const Teuchos::RCP< const panzer::ConnManager< LO, GO >> &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...
Teuchos::RCP< Tpetra::CrsMatrix< double, LO, GO, Kokkos::Compat::KokkosDeviceWrapperNode< PHX::Device > > > buildRHSMatrix(const panzer::UniqueGlobalIndexer< LO, GO > &sourceDOFManager, const Teuchos::RCP< const Tpetra::Map< LO, GO, Kokkos::Compat::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, one dimension of gradient (for hgrad basis), or one dimension of a vector field onto the target scalar basis. If you wish to project all values of a vector field or all the gradients of a scalar field, then you will need three separate RHS matrices to form the RHS for each independently. The vectors must be independent Tpetra vectors to solve multiple right hand sides with the linear solver.
const std::string & getType() const
Get type of basis.
Teuchos::RCP< Tpetra::CrsMatrix< double, LO, GO, Kokkos::Compat::KokkosDeviceWrapperNode< PHX::Device > > > buildMassMatrix()
Allocates, fills and returns a mass matrix for L2 projection onto a target basis. ...
Teuchos::RCP< CrsMatrixType > getGhostedTpetraMatrix(int i, int j) const
Kokkos::View< const int *, PHX::Device > offsets
Workset size is set to the total number of local elements in the MPI process.
Teuchos::RCP< CrsMatrixType > getTpetraMatrix(int i, int j) const