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>
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();
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>
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] = Kokkos::View<int*,PHX::Device>(
"offsets",offsets.size());
141 for(std::size_t i=0;i<offsets.size();i++)
142 scratch_offsets_[fd](i) = offsets[i];
145 scratch_lids_ = Kokkos::View<LO**,PHX::Device>(
"lids",gatherFields_[0].extent(0),
146 globalIndexer_->getElementBlockGIDCount(blockId));
148 indexerNames_.clear();
152 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
159 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(d.gedc->getDataObject(globalDataKey_));
161 if(tpetraContainer_==Teuchos::null) {
164 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(loc);
169 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
176 std::string blockId = this->wda(workset).block_id;
177 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
180 if (useTimeDerivativeSolutionVector_)
181 x = tpetraContainer_->get_dxdt();
183 x = tpetraContainer_->get_x();
185 auto x_data = x->template getLocalView<PHX::Device>();
187 globalIndexer_->getElementLIDs(this->wda(workset).cell_local_ids_k,scratch_lids_);
196 auto lids = scratch_lids_;
197 for (std::size_t fieldIndex=0; fieldIndex<gatherFields_.size();fieldIndex++) {
198 auto offsets = scratch_offsets_[fieldIndex];
199 auto gather_field = gatherFields_[fieldIndex];
201 Kokkos::parallel_for(localCellIds.size(), KOKKOS_LAMBDA (std::size_t worksetCellIndex) {
203 for(std::size_t basis=0;basis<
offsets.extent(0);basis++) {
205 LO lid =
lids(worksetCellIndex,offset);
208 gather_field(worksetCellIndex,basis) = x_data(lid,0);
218 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
223 : globalIndexer_(indexer)
224 , has_tangent_fields_(false)
226 typedef std::vector< std::vector<std::string> > vvstring;
231 const std::vector<std::string> & names = input.
getDofNames();
240 gatherFields_.resize(names.size());
241 for (std::size_t fd = 0; fd < names.size(); ++fd) {
243 PHX::MDField<ScalarT,Cell,NODE>(names[fd],basis->
functional);
244 this->addEvaluatedField(gatherFields_[fd]);
248 if (tangent_field_names.size()>0) {
249 TEUCHOS_ASSERT(gatherFields_.size() == tangent_field_names.size());
251 has_tangent_fields_ =
true;
252 tangentFields_.resize(gatherFields_.size());
253 for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
254 tangentFields_[fd].resize(tangent_field_names[fd].size());
255 for (std::size_t i=0; i<tangent_field_names[fd].size(); ++i) {
256 tangentFields_[fd][i] =
257 PHX::MDField<const RealT,Cell,NODE>(tangent_field_names[fd][i],basis->
functional);
258 this->addDependentField(tangentFields_[fd][i]);
264 std::string firstName =
"<none>";
266 firstName = names[0];
268 std::string n =
"GatherSolution (Tpetra): "+firstName+
" (Tangent)";
273 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
280 fieldIds_.resize(gatherFields_.size());
282 for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
283 const std::string& fieldName = indexerNames_[fd];
284 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
287 indexerNames_.clear();
291 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
298 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(d.gedc->getDataObject(globalDataKey_));
300 if(tpetraContainer_==Teuchos::null) {
303 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(loc);
308 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
315 std::string blockId = this->wda(workset).block_id;
316 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
319 if (useTimeDerivativeSolutionVector_)
320 x = tpetraContainer_->get_dxdt();
322 x = tpetraContainer_->get_x();
331 typedef typename PHX::MDField<ScalarT,Cell,NODE>::array_type::reference_type reference_type;
333 if (has_tangent_fields_) {
335 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
336 std::size_t cellLocalId = localCellIds[worksetCellIndex];
338 auto LIDs = globalIndexer_->getElementLIDs(cellLocalId);
341 for (std::size_t fieldIndex=0; fieldIndex<gatherFields_.size();fieldIndex++) {
342 int fieldNum = fieldIds_[fieldIndex];
343 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
345 const std::vector< PHX::MDField<const RealT,Cell,NODE> >& tf_ref =
346 tangentFields_[fieldIndex];
347 const std::size_t num_tf = tf_ref.size();
350 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
351 int offset = elmtOffset[basis];
352 LO lid = LIDs[offset];
353 reference_type gf_ref = (gatherFields_[fieldIndex])(worksetCellIndex,basis);
354 gf_ref.val() = x_array[lid];
355 for (std::size_t i=0; i<num_tf; ++i)
356 gf_ref.fastAccessDx(i) = tf_ref[i](worksetCellIndex,basis);
363 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
364 std::size_t cellLocalId = localCellIds[worksetCellIndex];
366 auto LIDs = globalIndexer_->getElementLIDs(cellLocalId);
369 for (std::size_t fieldIndex=0; fieldIndex<gatherFields_.size();fieldIndex++) {
370 int fieldNum = fieldIds_[fieldIndex];
371 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
374 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
375 int offset = elmtOffset[basis];
376 LO lid = LIDs[offset];
377 reference_type gf_ref = (gatherFields_[fieldIndex])(worksetCellIndex,basis);
378 gf_ref.val() = x_array[lid];
389 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
394 : globalIndexer_(indexer)
401 const std::vector<std::string> & names = input.
getDofNames();
413 gatherFields_.resize(names.size());
414 scratch_offsets_.resize(names.size());
415 for (std::size_t fd = 0; fd < names.size(); ++fd) {
416 PHX::MDField<ScalarT,Cell,NODE> f(names[fd],basis->
functional);
417 gatherFields_[fd] = f;
418 this->addEvaluatedField(gatherFields_[fd]);
422 std::string firstName =
"<none>";
424 firstName = names[0];
427 if(disableSensitivities_) {
428 std::string n =
"GatherSolution (Tpetra, No Sensitivities): "+firstName+
" (Jacobian)";
432 std::string n =
"GatherSolution (Tpetra): "+firstName+
" ("+PHX::typeAsString<EvalT>()+
") ";
438 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
445 fieldIds_.resize(gatherFields_.size());
447 const Workset & workset_0 = (*d.worksets_)[0];
448 std::string blockId = this->wda(workset_0).
block_id;
450 for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
452 const std::string& fieldName = indexerNames_[fd];
453 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
455 int fieldNum = fieldIds_[fd];
456 const std::vector<int> &
offsets = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
457 scratch_offsets_[fd] = Kokkos::View<int*,PHX::Device>(
"offsets",offsets.size());
458 for(std::size_t i=0;i<offsets.size();i++)
459 scratch_offsets_[fd](i) = offsets[i];
462 scratch_lids_ = Kokkos::View<LO**,PHX::Device>(
"lids",gatherFields_[0].extent(0),
463 globalIndexer_->getElementBlockGIDCount(blockId));
465 indexerNames_.clear();
469 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
475 using Teuchos::rcp_dynamic_cast;
482 if(!disableSensitivities_) {
483 if(d.first_sensitivities_name==sensitivitiesName_)
484 applySensitivities_ =
true;
486 applySensitivities_ =
false;
489 applySensitivities_ =
false;
496 std::string post = useTimeDerivativeSolutionVector_ ?
" - Xdot" :
" - X";
497 if(d.gedc->containsDataObject(globalDataKey_+post)) {
498 ged = d.gedc->getDataObject(globalDataKey_+post);
500 RCP<RO_GED> ro_ged = rcp_dynamic_cast<RO_GED>(ged,
true);
502 x_vector = ro_ged->getGhostedVector_Tpetra();
504 x_vector->template sync<PHX::Device>();
509 ged = d.gedc->getDataObject(globalDataKey_);
513 RCP<LOC> tpetraContainer = rcp_dynamic_cast<LOC>(ged);
516 if(loc_pair!=Teuchos::null) {
519 tpetraContainer = rcp_dynamic_cast<LOC>(loc);
522 if(tpetraContainer!=Teuchos::null) {
523 if (useTimeDerivativeSolutionVector_)
524 x_vector = tpetraContainer->get_dxdt();
526 x_vector = tpetraContainer->get_x();
528 x_vector->template sync<PHX::Device>();
536 RCP<RO_GED> ro_ged = rcp_dynamic_cast<RO_GED>(ged,
true);
538 x_vector = ro_ged->getGhostedVector_Tpetra();
541 x_vector->template sync<PHX::Device>();
545 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
550 std::string blockId = this->wda(workset).block_id;
552 double seed_value = 0.0;
553 if (useTimeDerivativeSolutionVector_) {
554 seed_value = workset.alpha;
556 else if (gatherSeedIndex_<0) {
557 seed_value = workset.beta;
559 else if(!useTimeDerivativeSolutionVector_) {
560 seed_value = workset.gather_seeds[gatherSeedIndex_];
569 if(!applySensitivities_)
573 bool use_seed =
true;
577 globalIndexer_->getElementLIDs(this->wda(workset).cell_local_ids_k,scratch_lids_);
582 functor_data.x_data = x_vector->template getLocalView<PHX::Device>();
583 functor_data.seed_value = seed_value;
584 functor_data.lids = scratch_lids_;
587 for(std::size_t fieldIndex=0;
588 fieldIndex<gatherFields_.size();fieldIndex++) {
591 functor_data.offsets = scratch_offsets_[fieldIndex];
592 functor_data.field = gatherFields_[fieldIndex];
595 Kokkos::parallel_for(workset.num_cells,*
this);
597 Kokkos::parallel_for(Kokkos::RangePolicy<PHX::Device,NoSeed>(0,workset.num_cells),*
this);
602 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
603 KOKKOS_INLINE_FUNCTION
608 for(std::size_t basis=0;basis<functor_data.offsets.extent(0);basis++) {
609 int offset = functor_data.offsets(basis);
610 LO lid = functor_data.lids(worksetCellIndex,offset);
613 functor_data.field(worksetCellIndex,basis).val() = functor_data.x_data(lid,0);
614 functor_data.field(worksetCellIndex,basis).fastAccessDx(offset) = functor_data.seed_value;
619 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
620 KOKKOS_INLINE_FUNCTION
625 for(std::size_t basis=0;basis<functor_data.offsets.extent(0);basis++) {
626 int offset = functor_data.offsets(basis);
627 LO lid = functor_data.lids(worksetCellIndex,offset);
630 functor_data.field(worksetCellIndex,basis).val() = functor_data.x_data(lid,0);
Kokkos::View< const LO **, PHX::Device > lids
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Gathers solution values from the Newton solution vector into the nodal fields of the field manager...
#define TEUCHOS_ASSERT(assertion_test)
Teuchos::RCP< PHX::DataLayout > functional
<Cell,Basis> or <Cell,Basis>
Kokkos::View< const int *, PHX::Device > offsets