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"
64template <
typename EvalT,
typename TRAITS,
typename S,
typename LO,
typename GO,
typename NodeT>
67 const Teuchos::RCP<const BlockedDOFManager> & indexer,
68 const Teuchos::ParameterList& p)
70 const std::vector<std::string>& names =
71 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"DOF Names"));
73 Teuchos::RCP<panzer::PureBasis> basis =
74 p.get< Teuchos::RCP<panzer::PureBasis> >(
"Basis");
76 for (std::size_t fd = 0; fd < names.size(); ++fd) {
77 PHX::MDField<ScalarT,Cell,NODE>
field = PHX::MDField<ScalarT,Cell,NODE>(names[fd],basis->functional);
78 this->addEvaluatedField(
field.fieldTag());
81 this->setName(
"Gather Solution");
88template <
typename TRAITS,
typename S,
typename LO,
typename GO,
typename NodeT>
91 const Teuchos::RCP<const BlockedDOFManager> & indexer,
92 const Teuchos::ParameterList& p)
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();
102 Teuchos::RCP<const panzer::PureBasis> basis = input.
getBasis();
110 gatherFields_.resize(names.size());
111 for (std::size_t fd = 0; fd < names.size(); ++fd) {
113 PHX::MDField<ScalarT,Cell,NODE>(names[fd],basis->functional);
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] =
127 PHX::MDField<const ScalarT,Cell,NODE>(tangent_field_names[fd][i],basis->functional);
128 this->addDependentField(tangentFields_[fd][i]);
134 std::string firstName =
"<none>";
136 firstName = names[0];
138 std::string n =
"GatherSolution (BlockedTpetra): "+firstName+
" (Residual)";
143template <
typename TRAITS,
typename S,
typename LO,
typename GO,
typename NodeT>
148 TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size());
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();
187template <
typename TRAITS,
typename S,
typename LO,
typename GO,
typename NodeT>
192 blockedContainer_ = Teuchos::rcp_dynamic_cast<const ContainerType>(d.gedc->getDataObject(globalDataKey_),
true);
196template <
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;
207 RCP<ProductVectorBase<ScalarT>> thyraBlockSolution;
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);
246template <
typename TRAITS,
typename S,
typename LO,
typename GO,
typename NodeT>
249 const Teuchos::RCP<const BlockedDOFManager> & indexer,
250 const Teuchos::ParameterList& p)
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();
260 Teuchos::RCP<const panzer::PureBasis> basis = input.
getBasis();
268 gatherFields_.resize(names.size());
269 for (std::size_t fd = 0; fd < names.size(); ++fd) {
271 PHX::MDField<ScalarT,Cell,NODE>(names[fd],basis->functional);
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] =
285 PHX::MDField<const ScalarT,Cell,NODE>(tangent_field_names[fd][i],basis->functional);
286 this->addDependentField(tangentFields_[fd][i]);
292 std::string firstName =
"<none>";
294 firstName = names[0];
296 std::string n =
"GatherSolution (BlockedTpetra): "+firstName+
" (Tangent)";
301template <
typename TRAITS,
typename S,
typename LO,
typename GO,
typename NodeT>
306 TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size());
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();
320template <
typename TRAITS,
typename S,
typename LO,
typename GO,
typename NodeT>
325 blockedContainer_ = Teuchos::rcp_dynamic_cast<const ContainerType>(d.gedc->getDataObject(globalDataKey_),
true);
329template <
typename TRAITS,
typename S,
typename LO,
typename GO,
typename NodeT>
334 using Teuchos::ArrayRCP;
335 using Teuchos::ptrFromRef;
336 using Teuchos::rcp_dynamic_cast;
339 using Thyra::SpmdVectorBase;
342 Teuchos::FancyOStream out(Teuchos::rcpFromRef(std::cout));
343 out.setShowProcRank(
true);
344 out.setOutputToRootOnly(-1);
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;
353 Teuchos::RCP<ProductVectorBase<double> > x;
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++) {
369 RCP<const MapType> x_map = blockedContainer_->getMapForBlock(GIDs[i].first);
371 LIDs[i] = x_map->getLocalElement(GIDs[i].second);
375 Teuchos::ArrayRCP<const double> local_x;
376 for (std::size_t fieldIndex=0; fieldIndex<gatherFields_.size();fieldIndex++) {
377 int fieldNum = fieldIds_[fieldIndex];
378 int indexerId = gidIndexer_->getFieldBlock(fieldNum);
381 RCP<SpmdVectorBase<double> > block_x = rcp_dynamic_cast<SpmdVectorBase<double> >(x->getNonconstVectorBlock(indexerId));
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();
408template <
typename TRAITS,
typename S,
typename LO,
typename GO,
typename NodeT>
411 const Teuchos::RCP<const BlockedDOFManager> & indexer,
412 const Teuchos::ParameterList& p)
413 : globalIndexer_(indexer)
418 const std::vector<std::string> & names = input.
getDofNames();
419 Teuchos::RCP<const panzer::PureBasis> basis = input.
getBasis();
427 gatherFields_.resize(names.size());
428 for (std::size_t fd = 0; fd < names.size(); ++fd) {
429 PHX::MDField<ScalarT,Cell,NODE> f(names[fd],basis->functional);
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>()+
") ";
451template <
typename TRAITS,
typename S,
typename LO,
typename GO,
typename NodeT>
456 TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size());
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)
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();
505template <
typename TRAITS,
typename S,
typename LO,
typename GO,
typename NodeT>
510 blockedContainer_ = Teuchos::rcp_dynamic_cast<const ContainerType>(d.gedc->getDataObject(globalDataKey_),
true);
514template <
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;
526 RCP<ProductVectorBase<double>> blockedSolution;
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;
PHX::View< const int * > offsets
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.
TRAITS::RealType RealType
Gathers solution values from the Newton solution vector into the nodal fields of the field manager.
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &vm)
void evaluateFields(typename TRAITS::EvalData d)
GatherSolution_BlockedTpetra(const Teuchos::RCP< const BlockedDOFManager > &indexer)
std::string block_id
DEPRECATED - use: getElementBlock()