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" 56 #include "Teuchos_FancyOStream.hpp" 58 #include "Tpetra_Vector.hpp" 59 #include "Tpetra_Map.hpp" 65 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
69 const Teuchos::ParameterList& p)
70 : globalIndexer_(indexer)
71 , has_tangent_fields_(false)
73 typedef std::vector< std::vector<std::string> > vvstring;
76 input.setParameterList(p);
78 const std::vector<std::string> & names =
input.getDofNames();
79 Teuchos::RCP<const panzer::PureBasis>
basis =
input.getBasis();
80 const vvstring & tangent_field_names =
input.getTangentNames();
82 indexerNames_ =
input.getIndexerNames();
83 useTimeDerivativeSolutionVector_ =
input.useTimeDerivativeSolutionVector();
84 globalDataKey_ =
input.getGlobalDataKey();
87 gatherFields_.resize(names.size());
88 for (std::size_t fd = 0; fd < names.size(); ++fd) {
90 PHX::MDField<ScalarT,Cell,NODE>(names[fd],
basis->functional);
91 this->addEvaluatedField(gatherFields_[fd]);
95 if (tangent_field_names.size()>0) {
96 TEUCHOS_ASSERT(gatherFields_.size() == tangent_field_names.size());
98 has_tangent_fields_ =
true;
99 tangentFields_.resize(gatherFields_.size());
100 for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
101 tangentFields_[fd].resize(tangent_field_names[fd].size());
102 for (std::size_t i=0; i<tangent_field_names[fd].size(); ++i) {
103 tangentFields_[fd][i] =
104 PHX::MDField<const ScalarT,Cell,NODE>(tangent_field_names[fd][i],
basis->functional);
105 this->addDependentField(tangentFields_[fd][i]);
111 std::string firstName =
"<none>";
113 firstName = names[0];
115 std::string n =
"GatherSolution (Tpetra): "+firstName+
" (Residual)";
120 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
125 TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size());
127 fieldIds_.resize(gatherFields_.size());
129 for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
130 const std::string& fieldName = indexerNames_[fd];
131 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
134 this->utils.setFieldData(gatherFields_[fd],fm);
137 if (has_tangent_fields_) {
138 for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd)
139 for (std::size_t i=0; i<tangentFields_[fd].size(); ++i)
140 this->utils.setFieldData(tangentFields_[fd][i],fm);
143 indexerNames_.clear();
147 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
154 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(d.gedc.getDataObject(globalDataKey_));
156 if(tpetraContainer_==Teuchos::null) {
158 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<
LOCPair_GlobalEvaluationData>(d.gedc.getDataObject(globalDataKey_),
true)->getGhostedLOC();
159 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(loc);
164 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
170 std::vector<LO> LIDs;
173 std::string blockId = this->wda(workset).block_id;
174 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
176 Teuchos::RCP<typename LOC::VectorType> x;
177 if (useTimeDerivativeSolutionVector_)
180 x = tpetraContainer_->get_x();
182 Teuchos::ArrayRCP<const double> x_array = x->get1dView();
190 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
191 std::size_t cellLocalId = localCellIds[worksetCellIndex];
193 LIDs = globalIndexer_->getElementLIDs(cellLocalId);
196 for (std::size_t fieldIndex=0; fieldIndex<gatherFields_.size();fieldIndex++) {
197 int fieldNum = fieldIds_[fieldIndex];
198 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
204 (gatherFields_[fieldIndex])(worksetCellIndex,
basis) = x_array[lid];
214 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
218 const Teuchos::ParameterList& p)
219 : globalIndexer_(indexer)
220 , has_tangent_fields_(false)
222 typedef std::vector< std::vector<std::string> > vvstring;
225 input.setParameterList(p);
227 const std::vector<std::string> & names =
input.getDofNames();
228 Teuchos::RCP<const panzer::PureBasis>
basis =
input.getBasis();
229 const vvstring & tangent_field_names =
input.getTangentNames();
231 indexerNames_ =
input.getIndexerNames();
232 useTimeDerivativeSolutionVector_ =
input.useTimeDerivativeSolutionVector();
233 globalDataKey_ =
input.getGlobalDataKey();
236 gatherFields_.resize(names.size());
237 for (std::size_t fd = 0; fd < names.size(); ++fd) {
239 PHX::MDField<ScalarT,Cell,NODE>(names[fd],
basis->functional);
240 this->addEvaluatedField(gatherFields_[fd]);
244 if (tangent_field_names.size()>0) {
245 TEUCHOS_ASSERT(gatherFields_.size() == tangent_field_names.size());
247 has_tangent_fields_ =
true;
248 tangentFields_.resize(gatherFields_.size());
249 for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
250 tangentFields_[fd].resize(tangent_field_names[fd].size());
251 for (std::size_t i=0; i<tangent_field_names[fd].size(); ++i) {
252 tangentFields_[fd][i] =
253 PHX::MDField<const RealT,Cell,NODE>(tangent_field_names[fd][i],
basis->functional);
254 this->addDependentField(tangentFields_[fd][i]);
260 std::string firstName =
"<none>";
262 firstName = names[0];
264 std::string n =
"GatherSolution (Tpetra): "+firstName+
" (Tangent)";
269 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
274 TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size());
276 fieldIds_.resize(gatherFields_.size());
278 for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
279 const std::string& fieldName = indexerNames_[fd];
280 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
283 this->utils.setFieldData(gatherFields_[fd],fm);
286 if (has_tangent_fields_) {
287 for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd)
288 for (std::size_t i=0; i<tangentFields_[fd].size(); ++i)
289 this->utils.setFieldData(tangentFields_[fd][i],fm);
292 indexerNames_.clear();
296 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
303 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(d.gedc.getDataObject(globalDataKey_));
305 if(tpetraContainer_==Teuchos::null) {
307 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<
LOCPair_GlobalEvaluationData>(d.gedc.getDataObject(globalDataKey_),
true)->getGhostedLOC();
308 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(loc);
313 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
319 std::vector<LO> LIDs;
322 std::string blockId = this->wda(workset).block_id;
323 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
325 Teuchos::RCP<typename LOC::VectorType> x;
326 if (useTimeDerivativeSolutionVector_)
329 x = tpetraContainer_->get_x();
331 Teuchos::ArrayRCP<const double> x_array = x->get1dView();
338 typedef typename PHX::MDField<ScalarT,Cell,NODE>::array_type::reference_type reference_type;
340 if (has_tangent_fields_) {
342 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
343 std::size_t cellLocalId = localCellIds[worksetCellIndex];
345 LIDs = globalIndexer_->getElementLIDs(cellLocalId);
348 for (std::size_t fieldIndex=0; fieldIndex<gatherFields_.size();fieldIndex++) {
349 int fieldNum = fieldIds_[fieldIndex];
350 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
352 const std::vector< PHX::MDField<const RealT,Cell,NODE> >& tf_ref =
353 tangentFields_[fieldIndex];
354 const std::size_t num_tf = tf_ref.size();
360 reference_type gf_ref = (gatherFields_[fieldIndex])(worksetCellIndex,
basis);
361 gf_ref.val() = x_array[lid];
362 for (std::size_t i=0; i<num_tf; ++i)
363 gf_ref.fastAccessDx(i) = tf_ref[i](worksetCellIndex,
basis);
370 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
371 std::size_t cellLocalId = localCellIds[worksetCellIndex];
373 LIDs = globalIndexer_->getElementLIDs(cellLocalId);
376 for (std::size_t fieldIndex=0; fieldIndex<gatherFields_.size();fieldIndex++) {
377 int fieldNum = fieldIds_[fieldIndex];
378 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
384 reference_type gf_ref = (gatherFields_[fieldIndex])(worksetCellIndex,
basis);
385 gf_ref.val() = x_array[lid];
396 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
400 const Teuchos::ParameterList& p)
401 : globalIndexer_(indexer)
406 input.setParameterList(p);
408 const std::vector<std::string> & names =
input.getDofNames();
409 Teuchos::RCP<const panzer::PureBasis>
basis =
input.getBasis();
412 indexerNames_ =
input.getIndexerNames();
413 useTimeDerivativeSolutionVector_ =
input.useTimeDerivativeSolutionVector();
414 globalDataKey_ =
input.getGlobalDataKey();
416 gatherSeedIndex_ =
input.getGatherSeedIndex();
417 sensitivitiesName_ =
input.getSensitivitiesName();
418 disableSensitivities_ = !
input.firstSensitivitiesAvailable();
420 gatherFields_.resize(names.size());
421 scratch_offsets_.resize(names.size());
422 for (std::size_t fd = 0; fd < names.size(); ++fd) {
423 PHX::MDField<ScalarT,Cell,NODE> f(names[fd],
basis->functional);
424 gatherFields_[fd] = f;
425 this->addEvaluatedField(gatherFields_[fd]);
429 std::string firstName =
"<none>";
431 firstName = names[0];
434 if(disableSensitivities_) {
435 std::string n =
"GatherSolution (Tpetra, No Sensitivities): "+firstName+
" (Jacobian)";
439 std::string n =
"GatherSolution (Tpetra): "+firstName+
" ("+PHX::typeAsString<EvalT>()+
") ";
445 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
450 TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size());
452 fieldIds_.resize(gatherFields_.size());
454 const Workset & workset_0 = (*d.worksets_)[0];
455 std::string blockId = this->wda(workset_0).
block_id;
457 for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
459 const std::string& fieldName = indexerNames_[fd];
460 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
463 this->utils.setFieldData(gatherFields_[fd],fm);
465 int fieldNum = fieldIds_[fd];
466 const std::vector<int> &
offsets = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
467 scratch_offsets_[fd] = Kokkos::View<int*,PHX::Device>(
"offsets",
offsets.size());
468 for(std::size_t i=0;i<
offsets.size();i++)
469 scratch_offsets_[fd](i) =
offsets[i];
472 scratch_lids_ = Kokkos::View<LO**,PHX::Device>(
"lids",gatherFields_[0].dimension_0(),
473 globalIndexer_->getElementBlockGIDCount(blockId));
475 indexerNames_.clear();
479 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
485 using Teuchos::rcp_dynamic_cast;
492 if(!disableSensitivities_) {
493 if(d.first_sensitivities_name==sensitivitiesName_)
494 applySensitivities_ =
true;
496 applySensitivities_ =
false;
499 applySensitivities_ =
false;
503 RCP<GlobalEvaluationData> ged;
506 std::string post = useTimeDerivativeSolutionVector_ ?
" - Xdot" :
" - X";
507 if(d.gedc.containsDataObject(globalDataKey_+post)) {
508 ged = d.gedc.getDataObject(globalDataKey_+post);
510 RCP<RO_GED> ro_ged = rcp_dynamic_cast<RO_GED>(ged,
true);
514 x_vector->template sync<PHX::Device>();
519 ged = d.gedc.getDataObject(globalDataKey_);
523 RCP<LOC> tpetraContainer = rcp_dynamic_cast<LOC>(ged);
526 if(loc_pair!=Teuchos::null) {
527 Teuchos::RCP<LinearObjContainer> loc = loc_pair->
getGhostedLOC();
529 tpetraContainer = rcp_dynamic_cast<LOC>(loc);
532 if(tpetraContainer!=Teuchos::null) {
533 if (useTimeDerivativeSolutionVector_)
534 x_vector = tpetraContainer->get_dxdt();
536 x_vector = tpetraContainer->get_x();
538 x_vector->template sync<PHX::Device>();
546 RCP<RO_GED> ro_ged = rcp_dynamic_cast<RO_GED>(ged,
true);
548 x_vector = ro_ged->getGhostedVector_Tpetra();
551 x_vector->template sync<PHX::Device>();
555 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
560 std::string blockId = this->wda(workset).block_id;
562 double seed_value = 0.0;
563 if (useTimeDerivativeSolutionVector_) {
564 seed_value = workset.alpha;
566 else if (gatherSeedIndex_<0) {
567 seed_value = workset.beta;
569 else if(!useTimeDerivativeSolutionVector_) {
570 seed_value = workset.gather_seeds[gatherSeedIndex_];
573 TEUCHOS_ASSERT(
false);
579 if(!applySensitivities_)
583 bool use_seed =
true;
587 globalIndexer_->getElementLIDs(this->wda(workset).cell_local_ids_k,scratch_lids_);
592 functor_data.x_data = x_vector->template getLocalView<PHX::Device>();
593 functor_data.seed_value = seed_value;
594 functor_data.lids = scratch_lids_;
597 for(std::size_t fieldIndex=0;
598 fieldIndex<gatherFields_.size();fieldIndex++) {
601 functor_data.offsets = scratch_offsets_[fieldIndex];
602 functor_data.field = gatherFields_[fieldIndex];
605 Kokkos::parallel_for(workset.num_cells,*
this);
607 Kokkos::parallel_for(Kokkos::RangePolicy<PHX::Device,NoSeed>(0,workset.num_cells),*
this);
612 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
613 KOKKOS_INLINE_FUNCTION
618 for(std::size_t
basis=0;
basis<functor_data.offsets.dimension_0();
basis++) {
620 LO lid = functor_data.lids(worksetCellIndex,
offset);
623 functor_data.field(worksetCellIndex,
basis).val() = functor_data.x_data(lid,0);
624 functor_data.field(worksetCellIndex,
basis).fastAccessDx(
offset) = functor_data.seed_value;
629 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
630 KOKKOS_INLINE_FUNCTION
635 for(std::size_t
basis=0;
basis<functor_data.offsets.dimension_0();
basis++) {
637 LO lid = functor_data.lids(worksetCellIndex,
offset);
640 functor_data.field(worksetCellIndex,
basis).val() = functor_data.x_data(lid,0);
PHX::MDField< const ScalarT > input
void operator()(const FieldMultTag< NUM_FIELD_MULT > &, const std::size_t &cell) const
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)
Teuchos::RCP< LinearObjContainer > getGhostedLOC() const
const Teuchos::RCP< VectorType > get_dxdt() const
Teuchos::RCP< const panzer::PureBasis > basis
Interpolates basis DOF values to IP DOF values.
Kokkos::View< const int *, PHX::Device > offsets