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>
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");
88 template <
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)";
143 template <
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();
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;
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 fieldGlobalIndexers_[fieldIndex]->getElementLIDs(localCellIds,worksetLIDs_);
219 currentWorksetLIDSubBlock = productVectorBlockIndex_[fieldIndex];
222 const auto& tpetraSolution = *((rcp_dynamic_cast<Thyra::TpetraVector<ScalarT,LO,GO,NodeT>>(thyraBlockSolution->getNonconstVectorBlock(productVectorBlockIndex_[fieldIndex]),
true))->getTpetraVector());
223 const auto& kokkosSolution = tpetraSolution.getLocalViewDevice(Tpetra::Access::ReadOnly);
226 const auto& fieldOffsets = fieldOffsets_[fieldIndex];
227 const auto& worksetLIDs = worksetLIDs_;
228 const auto& fieldValues = gatherFields_[fieldIndex];
230 Kokkos::parallel_for(Kokkos::RangePolicy<PHX::Device>(0,workset.num_cells), KOKKOS_LAMBDA (
const int& cell) {
231 for(
int basis=0; basis < static_cast<int>(fieldOffsets.size()); ++basis) {
232 const int lid = worksetLIDs(cell,fieldOffsets(basis));
233 fieldValues(cell,basis) = kokkosSolution(lid,0);
244 template <
typename TRAITS,
typename S,
typename LO,
typename GO,
typename NodeT>
247 const Teuchos::RCP<const BlockedDOFManager> & indexer,
248 const Teuchos::ParameterList& p)
249 : gidIndexer_(indexer)
250 , has_tangent_fields_(false)
252 typedef std::vector< std::vector<std::string> > vvstring;
257 const std::vector<std::string> & names = input.
getDofNames();
258 Teuchos::RCP<const panzer::PureBasis> basis = input.
getBasis();
266 gatherFields_.resize(names.size());
267 for (std::size_t fd = 0; fd < names.size(); ++fd) {
269 PHX::MDField<ScalarT,Cell,NODE>(names[fd],basis->functional);
270 this->addEvaluatedField(gatherFields_[fd]);
274 if (tangent_field_names.size()>0) {
275 TEUCHOS_ASSERT(gatherFields_.size() == tangent_field_names.size());
277 has_tangent_fields_ =
true;
278 tangentFields_.resize(gatherFields_.size());
279 for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
280 tangentFields_[fd].resize(tangent_field_names[fd].size());
281 for (std::size_t i=0; i<tangent_field_names[fd].size(); ++i) {
282 tangentFields_[fd][i] =
283 PHX::MDField<const ScalarT,Cell,NODE>(tangent_field_names[fd][i],basis->functional);
284 this->addDependentField(tangentFields_[fd][i]);
290 std::string firstName =
"<none>";
292 firstName = names[0];
294 std::string n =
"GatherSolution (BlockedTpetra): "+firstName+
" (Tangent)";
299 template <
typename TRAITS,
typename S,
typename LO,
typename GO,
typename NodeT>
304 TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size());
306 fieldIds_.resize(gatherFields_.size());
308 for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
310 const std::string& fieldName = indexerNames_[fd];
311 fieldIds_[fd] = gidIndexer_->getFieldNum(fieldName);
314 indexerNames_.clear();
318 template <
typename TRAITS,
typename S,
typename LO,
typename GO,
typename NodeT>
323 blockedContainer_ = Teuchos::rcp_dynamic_cast<
const ContainerType>(d.gedc->getDataObject(globalDataKey_),
true);
327 template <
typename TRAITS,
typename S,
typename LO,
typename GO,
typename NodeT>
332 using Teuchos::ArrayRCP;
333 using Teuchos::ptrFromRef;
334 using Teuchos::rcp_dynamic_cast;
337 using Thyra::SpmdVectorBase;
340 Teuchos::FancyOStream out(Teuchos::rcpFromRef(std::cout));
341 out.setShowProcRank(
true);
342 out.setOutputToRootOnly(-1);
344 std::vector<std::pair<int,GO> > GIDs;
345 std::vector<LO> LIDs;
348 std::string blockId = this->wda(workset).block_id;
349 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
351 Teuchos::RCP<ProductVectorBase<double> > x;
352 if (useTimeDerivativeSolutionVector_)
353 x = rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer_->get_dxdt());
355 x = rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer_->get_x());
358 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
359 LO cellLocalId = localCellIds[worksetCellIndex];
361 gidIndexer_->getElementGIDsPair(cellLocalId,GIDs,blockId);
364 LIDs.resize(GIDs.size());
365 for(std::size_t i=0;i<GIDs.size();i++) {
367 RCP<const MapType> x_map = blockedContainer_->getMapForBlock(GIDs[i].first);
369 LIDs[i] = x_map->getLocalElement(GIDs[i].second);
373 Teuchos::ArrayRCP<const double> local_x;
374 for (std::size_t fieldIndex=0; fieldIndex<gatherFields_.size();fieldIndex++) {
375 int fieldNum = fieldIds_[fieldIndex];
376 int indexerId = gidIndexer_->getFieldBlock(fieldNum);
379 RCP<SpmdVectorBase<double> > block_x = rcp_dynamic_cast<SpmdVectorBase<double> >(x->getNonconstVectorBlock(indexerId));
380 block_x->getLocalData(ptrFromRef(local_x));
382 const std::vector<int> & elmtOffset = gidIndexer_->getGIDFieldOffsets(blockId,fieldNum);
385 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
386 int offset = elmtOffset[basis];
387 int lid = LIDs[offset];
389 if (!has_tangent_fields_)
390 (gatherFields_[fieldIndex])(worksetCellIndex,basis) = local_x[lid];
392 (gatherFields_[fieldIndex])(worksetCellIndex,basis).val() = local_x[lid];
393 for (std::size_t i=0; i<tangentFields_[fieldIndex].size(); ++i)
394 (gatherFields_[fieldIndex])(worksetCellIndex,basis).fastAccessDx(i) =
395 tangentFields_[fieldIndex][i](worksetCellIndex,basis).val();
406 template <
typename TRAITS,
typename S,
typename LO,
typename GO,
typename NodeT>
409 const Teuchos::RCP<const BlockedDOFManager> & indexer,
410 const Teuchos::ParameterList& p)
411 : globalIndexer_(indexer)
416 const std::vector<std::string> & names = input.
getDofNames();
417 Teuchos::RCP<const panzer::PureBasis> basis = input.
getBasis();
425 gatherFields_.resize(names.size());
426 for (std::size_t fd = 0; fd < names.size(); ++fd) {
427 PHX::MDField<ScalarT,Cell,NODE> f(names[fd],basis->functional);
428 gatherFields_[fd] = f;
429 this->addEvaluatedField(gatherFields_[fd]);
433 std::string firstName =
"<none>";
435 firstName = names[0];
438 if(disableSensitivities_) {
439 std::string n =
"GatherSolution (BlockedTpetra, No Sensitivities): "+firstName+
" (Jacobian)";
443 std::string n =
"GatherSolution (BlockedTpetra): "+firstName+
" ("+PHX::print<EvalT>()+
") ";
449 template <
typename TRAITS,
typename S,
typename LO,
typename GO,
typename NodeT>
454 TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size());
456 const Workset & workset_0 = (*d.worksets_)[0];
457 const std::string blockId = this->wda(workset_0).
block_id;
459 fieldIds_.resize(gatherFields_.size());
460 fieldOffsets_.resize(gatherFields_.size());
461 productVectorBlockIndex_.resize(gatherFields_.size());
462 int maxElementBlockGIDCount = -1;
463 for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
465 const std::string fieldName = indexerNames_[fd];
466 const int globalFieldNum = globalIndexer_->getFieldNum(fieldName);
467 productVectorBlockIndex_[fd] = globalIndexer_->getFieldBlock(globalFieldNum);
468 const auto& subGlobalIndexer = globalIndexer_->getFieldDOFManagers()[productVectorBlockIndex_[fd]];
469 fieldIds_[fd] = subGlobalIndexer->getFieldNum(fieldName);
471 const std::vector<int>&
offsets = subGlobalIndexer->getGIDFieldOffsets(blockId,fieldIds_[fd]);
472 fieldOffsets_[fd] = PHX::View<int*>(
"GatherSolution_BlockedTpetra(Jacobian):fieldOffsets",
offsets.size());
473 auto hostOffsets = Kokkos::create_mirror_view(fieldOffsets_[fd]);
474 for (std::size_t i=0; i <
offsets.size(); ++i)
476 Kokkos::deep_copy(fieldOffsets_[fd], hostOffsets);
477 maxElementBlockGIDCount = std::max(subGlobalIndexer->getElementBlockGIDCount(blockId),maxElementBlockGIDCount);
483 worksetLIDs_ = PHX::View<LO**>(
"ScatterResidual_BlockedTpetra(Residual):worksetLIDs",
484 gatherFields_[0].extent(0),
485 maxElementBlockGIDCount);
488 const auto& blockGlobalIndexers = globalIndexer_->getFieldDOFManagers();
489 const int numBlocks =
static_cast<int>(globalIndexer_->getFieldDOFManagers().size());
490 blockOffsets_ = PHX::View<LO*>(
"GatherSolution_BlockedTpetra(Jacobian):blockOffsets_",
492 const auto hostBlockOffsets = Kokkos::create_mirror_view(blockOffsets_);
493 for (
int blk=0;blk<numBlocks;++blk) {
494 int blockOffset = globalIndexer_->getBlockGIDOffset(blockId,blk);
495 hostBlockOffsets(blk) = blockOffset;
497 hostBlockOffsets(numBlocks) = hostBlockOffsets(numBlocks-1) + blockGlobalIndexers[blockGlobalIndexers.size()-1]->getElementBlockGIDCount(blockId);
498 Kokkos::deep_copy(blockOffsets_,hostBlockOffsets);
500 indexerNames_.clear();
503 template <
typename TRAITS,
typename S,
typename LO,
typename GO,
typename NodeT>
508 blockedContainer_ = Teuchos::rcp_dynamic_cast<
const ContainerType>(d.gedc->getDataObject(globalDataKey_),
true);
512 template <
typename TRAITS,
typename S,
typename LO,
typename GO,
typename NodeT>
517 using Teuchos::rcp_dynamic_cast;
521 const auto& localCellIds = this->wda(workset).cell_local_ids_k;
524 RCP<ProductVectorBase<double>> blockedSolution;
525 if (useTimeDerivativeSolutionVector_) {
526 blockedSolution = rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer_->get_dxdt());
527 seedValue = workset.alpha;
530 blockedSolution = rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer_->get_x());
531 seedValue = workset.beta;
537 if(disableSensitivities_)
540 const int numFieldBlocks = globalIndexer_->getNumFieldBlocks();
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 blockIndexer->getElementLIDs(localCellIds,worksetLIDs_);
549 currentWorksetLIDSubBlock = productVectorBlockIndex_[fieldIndex];
552 const int blockRowIndex = productVectorBlockIndex_[fieldIndex];
553 const auto& subblockSolution = *((rcp_dynamic_cast<Thyra::TpetraVector<RealType,LO,GO,NodeT>>(blockedSolution->getNonconstVectorBlock(blockRowIndex),
true))->getTpetraVector());
554 const auto kokkosSolution = subblockSolution.getLocalViewDevice(Tpetra::Access::ReadOnly);
557 const PHX::View<const int*> fieldOffsets = fieldOffsets_[fieldIndex];
558 const PHX::View<const LO**> worksetLIDs = worksetLIDs_;
559 const PHX::View<ScalarT**> fieldValues = gatherFields_[fieldIndex].get_static_view();
560 const PHX::View<const LO*> blockOffsets = blockOffsets_;
561 auto blockOffsets_h = Kokkos::create_mirror_view(blockOffsets);
562 Kokkos::deep_copy(blockOffsets_h, blockOffsets);
563 const int blockStart = blockOffsets_h(blockRowIndex);
564 const int numDerivatives = blockOffsets_h(numFieldBlocks);
566 Kokkos::parallel_for(Kokkos::RangePolicy<PHX::Device>(0,workset.num_cells), KOKKOS_LAMBDA (
const int& cell) {
567 for (
int basis=0; basis < static_cast<int>(fieldOffsets.size()); ++basis) {
568 const int rowLID = worksetLIDs(cell,fieldOffsets(basis));
569 fieldValues(cell,basis) =
ScalarT(numDerivatives,kokkosSolution(rowLID,0));
570 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)
PHX::View< const int * > offsets
Gathers solution values from the Newton solution vector into the nodal fields of the field manager...
panzer::Traits::Jacobian::ScalarT ScalarT
TRAITS::RealType RealType
void evaluateFields(typename TRAITS::EvalData d)
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...