43 #ifndef PANZER_GATHER_SOLUTION_BLOCKED_EPETRA_IMPL_HPP
44 #define PANZER_GATHER_SOLUTION_BLOCKED_EPETRA_IMPL_HPP
46 #include "Teuchos_Assert.hpp"
47 #include "Phalanx_DataLayout.hpp"
52 #include "Panzer_TpetraLinearObjFactory.hpp"
57 #include "Teuchos_FancyOStream.hpp"
59 #include "Thyra_SpmdVectorBase.hpp"
60 #include "Thyra_ProductVectorBase.hpp"
62 #include "Tpetra_Map.hpp"
64 template <
typename EvalT,
typename TRAITS,
typename S,
typename LO,
typename GO,
typename NodeT>
70 const std::vector<std::string>& names =
76 for (std::size_t fd = 0; fd < names.size(); ++fd) {
78 this->addEvaluatedField(field.fieldTag());
81 this->setName(
"Gather Solution");
88 template <
typename TRAITS,
typename S,
typename LO,
typename GO,
typename NodeT>
93 : globalIndexer_(indexer)
94 , has_tangent_fields_(false)
96 typedef std::vector< std::vector<std::string> > vvstring;
101 const std::vector<std::string> & names = input.
getDofNames();
110 gatherFields_.resize(names.size());
111 for (std::size_t fd = 0; fd < names.size(); ++fd) {
114 this->addEvaluatedField(gatherFields_[fd]);
118 if (tangent_field_names.size()>0) {
119 TEUCHOS_ASSERT(gatherFields_.size() == tangent_field_names.size());
121 has_tangent_fields_ =
true;
122 tangentFields_.resize(gatherFields_.size());
123 for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
124 tangentFields_[fd].resize(tangent_field_names[fd].size());
125 for (std::size_t i=0; i<tangent_field_names[fd].size(); ++i) {
126 tangentFields_[fd][i] =
128 this->addDependentField(tangentFields_[fd][i]);
134 std::string firstName =
"<none>";
136 firstName = names[0];
138 std::string n =
"GatherSolution (BlockedTpetra): "+firstName+
" (Residual)";
143 template <
typename TRAITS,
typename S,
typename LO,
typename GO,
typename NodeT>
150 const Workset & workset_0 = (*d.worksets_)[0];
151 const std::string blockId = this->wda(workset_0).
block_id;
153 fieldIds_.resize(gatherFields_.size());
154 fieldOffsets_.resize(gatherFields_.size());
155 fieldGlobalIndexers_.resize(gatherFields_.size());
156 productVectorBlockIndex_.resize(gatherFields_.size());
157 int maxElementBlockGIDCount = -1;
158 for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
160 const std::string& fieldName = indexerNames_[fd];
161 const int globalFieldNum = globalIndexer_->getFieldNum(fieldName);
162 productVectorBlockIndex_[fd] = globalIndexer_->getFieldBlock(globalFieldNum);
163 fieldGlobalIndexers_[fd] = globalIndexer_->getFieldDOFManagers()[productVectorBlockIndex_[fd]];
164 fieldIds_[fd] = fieldGlobalIndexers_[fd]->getFieldNum(fieldName);
166 const std::vector<int>&
offsets = fieldGlobalIndexers_[fd]->getGIDFieldOffsets(blockId,fieldIds_[fd]);
167 fieldOffsets_[fd] = PHX::View<int*>(
"GatherSolution_BlockedTpetra(Residual):fieldOffsets",offsets.size());
168 auto hostFieldOffsets = Kokkos::create_mirror_view(fieldOffsets_[fd]);
169 for(std::size_t i=0; i < offsets.size(); ++i)
170 hostFieldOffsets(i) = offsets[i];
171 Kokkos::deep_copy(fieldOffsets_[fd],hostFieldOffsets);
173 maxElementBlockGIDCount = std::max(fieldGlobalIndexers_[fd]->getElementBlockGIDCount(blockId),maxElementBlockGIDCount);
179 worksetLIDs_ = PHX::View<LO**>(
"GatherSolution_BlockedTpetra(Residual):worksetLIDs",
180 gatherFields_[0].extent(0),
181 maxElementBlockGIDCount);
183 indexerNames_.clear();
187 template <
typename TRAITS,
typename S,
typename LO,
typename GO,
typename NodeT>
192 blockedContainer_ = Teuchos::rcp_dynamic_cast<
const ContainerType>(d.gedc->getDataObject(globalDataKey_),
true);
196 template <
typename TRAITS,
typename S,
typename LO,
typename GO,
typename NodeT>
201 using Teuchos::rcp_dynamic_cast;
205 const PHX::View<const int*>& localCellIds = this->wda(workset).cell_local_ids_k;
208 if (useTimeDerivativeSolutionVector_)
209 thyraBlockSolution = rcp_dynamic_cast<ProductVectorBase<ScalarT>>(blockedContainer_->get_dxdt(),
true);
211 thyraBlockSolution = rcp_dynamic_cast<ProductVectorBase<ScalarT>>(blockedContainer_->get_x(),
true);
214 int currentWorksetLIDSubBlock = -1;
215 for (std::size_t fieldIndex = 0; fieldIndex < gatherFields_.size(); fieldIndex++) {
217 if (productVectorBlockIndex_[fieldIndex] != currentWorksetLIDSubBlock) {
218 const std::string blockId = this->wda(workset).block_id;
219 const int num_dofs = fieldGlobalIndexers_[fieldIndex]->getElementBlockGIDCount(blockId);
220 fieldGlobalIndexers_[fieldIndex]->getElementLIDs(localCellIds,worksetLIDs_,num_dofs);
221 currentWorksetLIDSubBlock = productVectorBlockIndex_[fieldIndex];
224 const auto& tpetraSolution = *((rcp_dynamic_cast<Thyra::TpetraVector<ScalarT,LO,GO,NodeT>>(thyraBlockSolution->getNonconstVectorBlock(productVectorBlockIndex_[fieldIndex]),
true))->getTpetraVector());
225 const auto& kokkosSolution = tpetraSolution.getLocalViewDevice(Tpetra::Access::ReadOnly);
228 const auto& fieldOffsets = fieldOffsets_[fieldIndex];
229 const auto& worksetLIDs = worksetLIDs_;
230 const auto& fieldValues = gatherFields_[fieldIndex];
232 Kokkos::parallel_for(Kokkos::RangePolicy<PHX::Device>(0,workset.num_cells), KOKKOS_LAMBDA (
const int& cell) {
233 for(
int basis=0; basis < static_cast<int>(fieldOffsets.size()); ++basis) {
234 const int lid = worksetLIDs(cell,fieldOffsets(basis));
235 fieldValues(cell,basis) = kokkosSolution(lid,0);
246 template <
typename TRAITS,
typename S,
typename LO,
typename GO,
typename NodeT>
251 : gidIndexer_(indexer)
252 , has_tangent_fields_(false)
254 typedef std::vector< std::vector<std::string> > vvstring;
259 const std::vector<std::string> & names = input.
getDofNames();
268 gatherFields_.resize(names.size());
269 for (std::size_t fd = 0; fd < names.size(); ++fd) {
272 this->addEvaluatedField(gatherFields_[fd]);
276 if (tangent_field_names.size()>0) {
277 TEUCHOS_ASSERT(gatherFields_.size() == tangent_field_names.size());
279 has_tangent_fields_ =
true;
280 tangentFields_.resize(gatherFields_.size());
281 for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
282 tangentFields_[fd].resize(tangent_field_names[fd].size());
283 for (std::size_t i=0; i<tangent_field_names[fd].size(); ++i) {
284 tangentFields_[fd][i] =
286 this->addDependentField(tangentFields_[fd][i]);
292 std::string firstName =
"<none>";
294 firstName = names[0];
296 std::string n =
"GatherSolution (BlockedTpetra): "+firstName+
" (Tangent)";
301 template <
typename TRAITS,
typename S,
typename LO,
typename GO,
typename NodeT>
308 fieldIds_.resize(gatherFields_.size());
310 for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
312 const std::string& fieldName = indexerNames_[fd];
313 fieldIds_[fd] = gidIndexer_->getFieldNum(fieldName);
316 indexerNames_.clear();
320 template <
typename TRAITS,
typename S,
typename LO,
typename GO,
typename NodeT>
325 blockedContainer_ = Teuchos::rcp_dynamic_cast<
const ContainerType>(d.gedc->getDataObject(globalDataKey_),
true);
329 template <
typename TRAITS,
typename S,
typename LO,
typename GO,
typename NodeT>
335 using Teuchos::ptrFromRef;
336 using Teuchos::rcp_dynamic_cast;
339 using Thyra::SpmdVectorBase;
346 std::vector<std::pair<int,GO> > GIDs;
347 std::vector<LO> LIDs;
350 std::string blockId = this->wda(workset).block_id;
351 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
354 if (useTimeDerivativeSolutionVector_)
355 x = rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer_->get_dxdt());
357 x = rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer_->get_x());
360 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
361 LO cellLocalId = localCellIds[worksetCellIndex];
363 gidIndexer_->getElementGIDsPair(cellLocalId,GIDs,blockId);
366 LIDs.resize(GIDs.size());
367 for(std::size_t i=0;i<GIDs.size();i++) {
371 LIDs[i] = x_map->getLocalElement(GIDs[i].second);
376 for (std::size_t fieldIndex=0; fieldIndex<gatherFields_.size();fieldIndex++) {
377 int fieldNum = fieldIds_[fieldIndex];
378 int indexerId = gidIndexer_->getFieldBlock(fieldNum);
382 block_x->getLocalData(ptrFromRef(local_x));
384 const std::vector<int> & elmtOffset = gidIndexer_->getGIDFieldOffsets(blockId,fieldNum);
387 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
388 int offset = elmtOffset[basis];
389 int lid = LIDs[offset];
391 if (!has_tangent_fields_)
392 (gatherFields_[fieldIndex])(worksetCellIndex,basis) = local_x[lid];
394 (gatherFields_[fieldIndex])(worksetCellIndex,basis).val() = local_x[lid];
395 for (std::size_t i=0; i<tangentFields_[fieldIndex].size(); ++i)
396 (gatherFields_[fieldIndex])(worksetCellIndex,basis).fastAccessDx(i) =
397 tangentFields_[fieldIndex][i](worksetCellIndex,basis).val();
408 template <
typename TRAITS,
typename S,
typename LO,
typename GO,
typename NodeT>
413 : globalIndexer_(indexer)
418 const std::vector<std::string> & names = input.
getDofNames();
427 gatherFields_.resize(names.size());
428 for (std::size_t fd = 0; fd < names.size(); ++fd) {
430 gatherFields_[fd] = f;
431 this->addEvaluatedField(gatherFields_[fd]);
435 std::string firstName =
"<none>";
437 firstName = names[0];
440 if(disableSensitivities_) {
441 std::string n =
"GatherSolution (BlockedTpetra, No Sensitivities): "+firstName+
" (Jacobian)";
445 std::string n =
"GatherSolution (BlockedTpetra): "+firstName+
" ("+PHX::print<EvalT>()+
") ";
451 template <
typename TRAITS,
typename S,
typename LO,
typename GO,
typename NodeT>
458 const Workset & workset_0 = (*d.worksets_)[0];
459 const std::string blockId = this->wda(workset_0).
block_id;
461 fieldIds_.resize(gatherFields_.size());
462 fieldOffsets_.resize(gatherFields_.size());
463 productVectorBlockIndex_.resize(gatherFields_.size());
464 int maxElementBlockGIDCount = -1;
465 for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
467 const std::string fieldName = indexerNames_[fd];
468 const int globalFieldNum = globalIndexer_->getFieldNum(fieldName);
469 productVectorBlockIndex_[fd] = globalIndexer_->getFieldBlock(globalFieldNum);
470 const auto& subGlobalIndexer = globalIndexer_->getFieldDOFManagers()[productVectorBlockIndex_[fd]];
471 fieldIds_[fd] = subGlobalIndexer->getFieldNum(fieldName);
473 const std::vector<int>&
offsets = subGlobalIndexer->getGIDFieldOffsets(blockId,fieldIds_[fd]);
474 fieldOffsets_[fd] = PHX::View<int*>(
"GatherSolution_BlockedTpetra(Jacobian):fieldOffsets",offsets.size());
475 auto hostOffsets = Kokkos::create_mirror_view(fieldOffsets_[fd]);
476 for (std::size_t i=0; i < offsets.size(); ++i)
477 hostOffsets(i) = offsets[i];
478 Kokkos::deep_copy(fieldOffsets_[fd], hostOffsets);
479 maxElementBlockGIDCount = std::max(subGlobalIndexer->getElementBlockGIDCount(blockId),maxElementBlockGIDCount);
485 worksetLIDs_ = PHX::View<LO**>(
"ScatterResidual_BlockedTpetra(Residual):worksetLIDs",
486 gatherFields_[0].extent(0),
487 maxElementBlockGIDCount);
490 const auto& blockGlobalIndexers = globalIndexer_->getFieldDOFManagers();
491 const int numBlocks =
static_cast<int>(globalIndexer_->getFieldDOFManagers().size());
492 blockOffsets_ = PHX::View<LO*>(
"GatherSolution_BlockedTpetra(Jacobian):blockOffsets_",
494 const auto hostBlockOffsets = Kokkos::create_mirror_view(blockOffsets_);
495 for (
int blk=0;blk<numBlocks;++blk) {
496 int blockOffset = globalIndexer_->getBlockGIDOffset(blockId,blk);
497 hostBlockOffsets(blk) = blockOffset;
499 hostBlockOffsets(numBlocks) = hostBlockOffsets(numBlocks-1) + blockGlobalIndexers[blockGlobalIndexers.size()-1]->getElementBlockGIDCount(blockId);
500 Kokkos::deep_copy(blockOffsets_,hostBlockOffsets);
502 indexerNames_.clear();
505 template <
typename TRAITS,
typename S,
typename LO,
typename GO,
typename NodeT>
510 blockedContainer_ = Teuchos::rcp_dynamic_cast<
const ContainerType>(d.gedc->getDataObject(globalDataKey_),
true);
514 template <
typename TRAITS,
typename S,
typename LO,
typename GO,
typename NodeT>
519 using Teuchos::rcp_dynamic_cast;
523 const auto& localCellIds = this->wda(workset).cell_local_ids_k;
527 if (useTimeDerivativeSolutionVector_) {
528 blockedSolution = rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer_->get_dxdt());
529 seedValue = workset.alpha;
532 blockedSolution = rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer_->get_x());
533 seedValue = workset.beta;
539 if(disableSensitivities_)
543 int currentWorksetLIDSubBlock = -1;
544 for (std::size_t fieldIndex = 0; fieldIndex < gatherFields_.size(); fieldIndex++) {
546 if (productVectorBlockIndex_[fieldIndex] != currentWorksetLIDSubBlock) {
547 const auto& blockIndexer = globalIndexer_->getFieldDOFManagers()[productVectorBlockIndex_[fieldIndex]];
548 const std::string blockId = this->wda(workset).block_id;
549 const int num_dofs = globalIndexer_->getFieldDOFManagers()[productVectorBlockIndex_[fieldIndex]]->getElementBlockGIDCount(blockId);
550 blockIndexer->getElementLIDs(localCellIds,worksetLIDs_,num_dofs);
551 currentWorksetLIDSubBlock = productVectorBlockIndex_[fieldIndex];
554 const int blockRowIndex = productVectorBlockIndex_[fieldIndex];
555 const auto& subblockSolution = *((rcp_dynamic_cast<Thyra::TpetraVector<RealType,LO,GO,NodeT>>(blockedSolution->getNonconstVectorBlock(blockRowIndex),
true))->getTpetraVector());
556 const auto kokkosSolution = subblockSolution.getLocalViewDevice(Tpetra::Access::ReadOnly);
559 const PHX::View<const int*> fieldOffsets = fieldOffsets_[fieldIndex];
560 const PHX::View<const LO**> worksetLIDs = worksetLIDs_;
561 const PHX::View<ScalarT**> fieldValues = gatherFields_[fieldIndex].get_static_view();
562 const PHX::View<const LO*> blockOffsets = blockOffsets_;
563 auto blockOffsets_h = Kokkos::create_mirror_view(blockOffsets);
564 Kokkos::deep_copy(blockOffsets_h, blockOffsets);
565 const int blockStart = blockOffsets_h(blockRowIndex);
567 Kokkos::parallel_for(Kokkos::RangePolicy<PHX::Device>(0,workset.num_cells), KOKKOS_LAMBDA (
const int& cell) {
568 for (
int basis=0; basis < static_cast<int>(fieldOffsets.size()); ++basis) {
569 const int rowLID = worksetLIDs(cell,fieldOffsets(basis));
570 fieldValues(cell,basis).zero();
571 fieldValues(cell,basis).val() = kokkosSolution(rowLID,0);
572 fieldValues(cell,basis).fastAccessDx(blockStart+fieldOffsets(basis)) = seedValue;
GatherSolution_BlockedTpetra(const Teuchos::RCP< const BlockedDOFManager > &indexer)
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &vm)
basic_FancyOStream & setShowProcRank(const bool showProcRank)
T & get(const std::string &name, T def_value)
PHX::View< const int * > offsets
TRAITS::RealType RealType
void evaluateFields(typename TRAITS::EvalData d)
basic_FancyOStream & setOutputToRootOnly(const int rootRank)
virtual void preEvaluate(typename Traits::PreEvalData d)=0
std::string block_id
DEPRECATED - use: getElementBlock()
PHX::MDField< ScalarT, panzer::Cell, panzer::BASIS > field
A field to which we'll contribute, or in which we'll store, the result of computing this integral...
#define TEUCHOS_ASSERT(assertion_test)
Teuchos::RCP< PHX::DataLayout > functional
<Cell,Basis> or <Cell,Basis>