43 #ifndef PANZER_GATHER_SOLUTION_TPETRA_IMPL_HPP 44 #define PANZER_GATHER_SOLUTION_TPETRA_IMPL_HPP 46 #include "Teuchos_Assert.hpp" 47 #include "Phalanx_DataLayout.hpp" 57 #include "Teuchos_FancyOStream.hpp" 59 #include "Tpetra_Vector.hpp" 60 #include "Tpetra_Map.hpp" 66 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
69 const Teuchos::RCP<const panzer::GlobalIndexer> & indexer,
70 const Teuchos::ParameterList& p)
71 : globalIndexer_(indexer)
72 , has_tangent_fields_(false)
74 typedef std::vector< std::vector<std::string> > vvstring;
79 const std::vector<std::string> & names = input.
getDofNames();
80 Teuchos::RCP<const panzer::PureBasis> basis = input.
getBasis();
88 gatherFields_.resize(names.size());
89 for (std::size_t fd = 0; fd < names.size(); ++fd) {
91 PHX::MDField<ScalarT,Cell,NODE>(names[fd],basis->functional);
92 this->addEvaluatedField(gatherFields_[fd]);
96 if (tangent_field_names.size()>0) {
97 TEUCHOS_ASSERT(gatherFields_.size() == tangent_field_names.size());
99 has_tangent_fields_ =
true;
100 tangentFields_.resize(gatherFields_.size());
101 for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
102 tangentFields_[fd].resize(tangent_field_names[fd].size());
103 for (std::size_t i=0; i<tangent_field_names[fd].size(); ++i) {
104 tangentFields_[fd][i] =
105 PHX::MDField<const ScalarT,Cell,NODE>(tangent_field_names[fd][i],basis->functional);
106 this->addDependentField(tangentFields_[fd][i]);
112 std::string firstName =
"<none>";
114 firstName = names[0];
116 std::string n =
"GatherSolution (Tpetra): "+firstName+
" (Residual)";
121 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
126 TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size());
128 fieldIds_.resize(gatherFields_.size());
130 const Workset & workset_0 = (*d.worksets_)[0];
131 std::string blockId = this->wda(workset_0).
block_id;
132 scratch_offsets_.resize(gatherFields_.size());
134 for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
135 const std::string& fieldName = indexerNames_[fd];
136 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
138 int fieldNum = fieldIds_[fd];
139 const std::vector<int> &
offsets = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
140 scratch_offsets_[fd] = PHX::View<int*>(
"offsets",
offsets.size());
141 Kokkos::deep_copy(scratch_offsets_[fd], Kokkos::View<const int*, Kokkos::HostSpace, Kokkos::MemoryUnmanaged>(
offsets.data(),
offsets.size()));
144 scratch_lids_ = PHX::View<LO**>(
"lids",gatherFields_[0].extent(0),
145 globalIndexer_->getElementBlockGIDCount(blockId));
147 indexerNames_.clear();
151 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
158 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(d.gedc->getDataObject(globalDataKey_));
160 if(tpetraContainer_==Teuchos::null) {
162 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<
LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),
true)->getGhostedLOC();
163 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(loc);
168 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
175 std::string blockId = this->wda(workset).block_id;
176 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
178 Teuchos::RCP<typename LOC::VectorType> x;
179 if (useTimeDerivativeSolutionVector_)
182 x = tpetraContainer_->get_x();
184 auto x_data = x->getLocalViewDevice(Tpetra::Access::ReadOnly);
186 globalIndexer_->getElementLIDs(this->wda(workset).cell_local_ids_k,scratch_lids_);
195 auto lids = scratch_lids_;
196 for (std::size_t fieldIndex=0; fieldIndex<gatherFields_.size();fieldIndex++) {
197 auto offsets = scratch_offsets_[fieldIndex];
198 auto gather_field = gatherFields_[fieldIndex];
200 Kokkos::parallel_for(localCellIds.size(), KOKKOS_LAMBDA (std::size_t worksetCellIndex) {
202 for(std::size_t basis=0;basis<
offsets.extent(0);basis++) {
204 LO lid =
lids(worksetCellIndex,offset);
207 gather_field(worksetCellIndex,basis) = x_data(lid,0);
217 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
220 const Teuchos::RCP<const panzer::GlobalIndexer> & indexer,
221 const Teuchos::ParameterList& p)
222 : globalIndexer_(indexer)
223 , has_tangent_fields_(false)
225 typedef std::vector< std::vector<std::string> > vvstring;
230 const std::vector<std::string> & names = input.
getDofNames();
231 Teuchos::RCP<const panzer::PureBasis> basis = input.
getBasis();
239 gatherFields_.resize(names.size());
240 for (std::size_t fd = 0; fd < names.size(); ++fd) {
242 PHX::MDField<ScalarT,Cell,NODE>(names[fd],basis->functional);
243 this->addEvaluatedField(gatherFields_[fd]);
246 this->addUnsharedField(gatherFields_[fd].fieldTag().clone());
250 if (tangent_field_names.size()>0) {
251 TEUCHOS_ASSERT(gatherFields_.size() == tangent_field_names.size());
253 has_tangent_fields_ =
true;
254 tangentFields_.resize(gatherFields_.size());
255 for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
256 tangentFields_[fd].resize(tangent_field_names[fd].size());
257 for (std::size_t i=0; i<tangent_field_names[fd].size(); ++i) {
258 tangentFields_[fd][i] =
259 PHX::MDField<const RealT,Cell,NODE>(tangent_field_names[fd][i],basis->functional);
260 this->addDependentField(tangentFields_[fd][i]);
266 std::string firstName =
"<none>";
268 firstName = names[0];
270 std::string n =
"GatherSolution (Tpetra): "+firstName+
" (Tangent)";
275 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
280 TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size());
282 fieldIds_.resize(gatherFields_.size());
284 for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
285 const std::string& fieldName = indexerNames_[fd];
286 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
289 indexerNames_.clear();
293 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
300 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(d.gedc->getDataObject(globalDataKey_));
302 if(tpetraContainer_==Teuchos::null) {
304 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<
LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),
true)->getGhostedLOC();
305 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(loc);
310 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
317 std::string blockId = this->wda(workset).block_id;
318 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
320 Teuchos::RCP<typename LOC::VectorType> x;
321 if (useTimeDerivativeSolutionVector_)
324 x = tpetraContainer_->get_x();
326 Teuchos::ArrayRCP<const double> x_array = x->get1dView();
333 typedef typename PHX::MDField<ScalarT,Cell,NODE>::array_type::reference_type reference_type;
335 if (has_tangent_fields_) {
337 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
338 std::size_t cellLocalId = localCellIds[worksetCellIndex];
340 auto LIDs = globalIndexer_->getElementLIDs(cellLocalId);
343 for (std::size_t fieldIndex=0; fieldIndex<gatherFields_.size();fieldIndex++) {
344 int fieldNum = fieldIds_[fieldIndex];
345 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
347 const std::vector< PHX::MDField<const RealT,Cell,NODE> >& tf_ref =
348 tangentFields_[fieldIndex];
349 const std::size_t num_tf = tf_ref.size();
352 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
353 int offset = elmtOffset[basis];
354 LO lid = LIDs[offset];
355 reference_type gf_ref = (gatherFields_[fieldIndex])(worksetCellIndex,basis);
356 gf_ref.val() = x_array[lid];
357 for (std::size_t i=0; i<num_tf; ++i)
358 gf_ref.fastAccessDx(i) = tf_ref[i](worksetCellIndex,basis);
365 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
366 std::size_t cellLocalId = localCellIds[worksetCellIndex];
368 auto LIDs = globalIndexer_->getElementLIDs(cellLocalId);
371 for (std::size_t fieldIndex=0; fieldIndex<gatherFields_.size();fieldIndex++) {
372 int fieldNum = fieldIds_[fieldIndex];
373 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
376 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
377 int offset = elmtOffset[basis];
378 LO lid = LIDs[offset];
379 reference_type gf_ref = (gatherFields_[fieldIndex])(worksetCellIndex,basis);
380 gf_ref.val() = x_array[lid];
391 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
394 const Teuchos::RCP<const panzer::GlobalIndexer> & indexer,
395 const Teuchos::ParameterList& p)
396 : globalIndexer_(indexer)
403 const std::vector<std::string> & names = input.
getDofNames();
404 Teuchos::RCP<const panzer::PureBasis> basis = input.
getBasis();
415 gatherFields_.resize(names.size());
416 scratch_offsets_.resize(names.size());
417 for (std::size_t fd = 0; fd < names.size(); ++fd) {
418 PHX::MDField<ScalarT,Cell,NODE> f(names[fd],basis->functional);
419 gatherFields_[fd] = f;
420 this->addEvaluatedField(gatherFields_[fd]);
423 this->addUnsharedField(gatherFields_[fd].fieldTag().clone());
427 std::string firstName =
"<none>";
429 firstName = names[0];
432 if(disableSensitivities_) {
433 std::string n =
"GatherSolution (Tpetra, No Sensitivities): "+firstName+
" (Jacobian)";
437 std::string n =
"GatherSolution (Tpetra): "+firstName+
" (Jacobian) ";
443 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
448 TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size());
450 fieldIds_.resize(gatherFields_.size());
452 const Workset & workset_0 = (*d.worksets_)[0];
453 std::string blockId = this->wda(workset_0).
block_id;
455 for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
457 const std::string& fieldName = indexerNames_[fd];
458 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
460 int fieldNum = fieldIds_[fd];
461 const std::vector<int> &
offsets = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
462 scratch_offsets_[fd] = PHX::View<int*>(
"offsets",
offsets.size());
463 Kokkos::deep_copy(scratch_offsets_[fd], Kokkos::View<const int*, Kokkos::HostSpace, Kokkos::MemoryUnmanaged>(
offsets.data(),
offsets.size()));
466 scratch_lids_ = PHX::View<LO**>(
"lids",gatherFields_[0].extent(0),
467 globalIndexer_->getElementBlockGIDCount(blockId));
469 indexerNames_.clear();
473 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
479 using Teuchos::rcp_dynamic_cast;
486 if(!disableSensitivities_) {
487 if(d.first_sensitivities_name==sensitivitiesName_)
488 applySensitivities_ =
true;
490 applySensitivities_ =
false;
493 applySensitivities_ =
false;
497 RCP<GlobalEvaluationData> ged;
500 std::string post = useTimeDerivativeSolutionVector_ ?
" - Xdot" :
" - X";
501 if(d.gedc->containsDataObject(globalDataKey_+post)) {
502 ged = d.gedc->getDataObject(globalDataKey_+post);
504 RCP<RO_GED> ro_ged = rcp_dynamic_cast<RO_GED>(ged,
true);
511 ged = d.gedc->getDataObject(globalDataKey_);
515 RCP<LOC> tpetraContainer = rcp_dynamic_cast<LOC>(ged);
518 if(loc_pair!=Teuchos::null) {
519 Teuchos::RCP<LinearObjContainer> loc = loc_pair->
getGhostedLOC();
521 tpetraContainer = rcp_dynamic_cast<LOC>(loc);
524 if(tpetraContainer!=Teuchos::null) {
525 if (useTimeDerivativeSolutionVector_)
526 x_vector = tpetraContainer->get_dxdt();
528 x_vector = tpetraContainer->get_x();
536 RCP<RO_GED> ro_ged = rcp_dynamic_cast<RO_GED>(ged,
true);
538 x_vector = ro_ged->getGhostedVector_Tpetra();
543 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
548 std::string blockId = this->wda(workset).block_id;
550 double seed_value = 0.0;
551 if (useTimeDerivativeSolutionVector_) {
552 seed_value = workset.alpha;
554 else if (gatherSeedIndex_<0) {
555 seed_value = workset.beta;
557 else if(!useTimeDerivativeSolutionVector_) {
558 seed_value = workset.gather_seeds[gatherSeedIndex_];
561 TEUCHOS_ASSERT(
false);
567 if(!applySensitivities_)
573 functor_data.dos = 0;
574 if (this->wda.getDetailsIndex() == 1)
577 functor_data.dos = globalIndexer_->getElementBlockGIDCount(workset.details(0).block_id);
581 bool use_seed =
true;
585 globalIndexer_->getElementLIDs(this->wda(workset).cell_local_ids_k,scratch_lids_);
590 functor_data.x_data = x_vector->getLocalViewDevice(Tpetra::Access::ReadOnly);
591 functor_data.seed_value = seed_value;
592 functor_data.lids = scratch_lids_;
595 for(std::size_t fieldIndex=0;
596 fieldIndex<gatherFields_.size();fieldIndex++) {
599 functor_data.offsets = scratch_offsets_[fieldIndex];
600 functor_data.field = gatherFields_[fieldIndex];
603 Kokkos::parallel_for(workset.num_cells,*
this);
605 Kokkos::parallel_for(Kokkos::RangePolicy<PHX::Device,NoSeed>(0,workset.num_cells),*
this);
607 functor_data.x_data = Kokkos::View<const double**, Kokkos::LayoutLeft,PHX::Device>();
611 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
612 KOKKOS_INLINE_FUNCTION
617 for(std::size_t basis=0;basis<functor_data.offsets.extent(0);basis++) {
618 int offset = functor_data.offsets(basis);
619 LO lid = functor_data.lids(worksetCellIndex,offset);
622 if (functor_data.dos == 0)
623 functor_data.field(worksetCellIndex,basis).val() = functor_data.x_data(lid,0);
625 functor_data.field(worksetCellIndex,basis) =
ScalarT(functor_data.x_data(lid,0));
627 functor_data.field(worksetCellIndex,basis).fastAccessDx(functor_data.dos + offset) = functor_data.seed_value;
632 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
633 KOKKOS_INLINE_FUNCTION
638 for(std::size_t basis=0;basis<functor_data.offsets.extent(0);basis++) {
639 int offset = functor_data.offsets(basis);
640 LO lid = functor_data.lids(worksetCellIndex,offset);
643 functor_data.field(worksetCellIndex,basis).val() = functor_data.x_data(lid,0);
panzer::Traits::Jacobian::ScalarT ScalarT
PHX::View< const int * > offsets
Gathers solution values from the Newton solution vector into the nodal fields of the field manager...
Teuchos::RCP< VectorType > getGhostedVector_Tpetra() const
Get the ghosted vector (Tpetra version)
std::string block_id
DEPRECATED - use: getElementBlock()
Teuchos::RCP< LinearObjContainer > getGhostedLOC() const
const Teuchos::RCP< VectorType > get_dxdt() const
PHX::View< const LO ** > lids