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"
58 #include "Teuchos_FancyOStream.hpp"
60 #include "Tpetra_Vector.hpp"
61 #include "Tpetra_Map.hpp"
67 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
72 : globalIndexer_(indexer)
73 , has_tangent_fields_(false)
75 typedef std::vector< std::vector<std::string> > vvstring;
80 const std::vector<std::string> & names = input.
getDofNames();
89 gatherFields_.resize(names.size());
90 for (std::size_t fd = 0; fd < names.size(); ++fd) {
93 this->addEvaluatedField(gatherFields_[fd]);
97 if (tangent_field_names.size()>0) {
98 TEUCHOS_ASSERT(gatherFields_.size() == tangent_field_names.size());
100 has_tangent_fields_ =
true;
101 tangentFields_.resize(gatherFields_.size());
102 for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
103 tangentFields_[fd].resize(tangent_field_names[fd].size());
104 for (std::size_t i=0; i<tangent_field_names[fd].size(); ++i) {
105 tangentFields_[fd][i] =
107 this->addDependentField(tangentFields_[fd][i]);
113 std::string firstName =
"<none>";
115 firstName = names[0];
117 std::string n =
"GatherSolution (Tpetra): "+firstName+
" (Residual)";
122 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
129 fieldIds_.resize(gatherFields_.size());
131 const Workset & workset_0 = (*d.worksets_)[0];
132 std::string blockId = this->wda(workset_0).
block_id;
133 scratch_offsets_.resize(gatherFields_.size());
135 for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
136 const std::string& fieldName = indexerNames_[fd];
137 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
139 int fieldNum = fieldIds_[fd];
140 const std::vector<int> &
offsets = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
141 scratch_offsets_[fd] = PHX::View<int*>(
"offsets",offsets.size());
145 scratch_lids_ = PHX::View<LO**>(
"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->getLocalViewDevice(Tpetra::Access::ReadOnly);
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].get_static_view();
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) {
244 this->addEvaluatedField(gatherFields_[fd]);
247 this->addUnsharedField(gatherFields_[fd].fieldTag().clone());
251 if (tangent_field_names.size()>0) {
252 TEUCHOS_ASSERT(gatherFields_.size() == tangent_field_names.size());
254 has_tangent_fields_ =
true;
255 tangentFields_.resize(gatherFields_.size());
256 for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
257 tangentFields_[fd].resize(tangent_field_names[fd].size());
258 for (std::size_t i=0; i<tangent_field_names[fd].size(); ++i) {
259 tangentFields_[fd][i] =
261 this->addDependentField(tangentFields_[fd][i]);
267 std::string firstName =
"<none>";
269 firstName = names[0];
271 std::string n =
"GatherSolution (Tpetra): "+firstName+
" (Tangent)";
276 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
283 fieldIds_.resize(gatherFields_.size());
289 tangentInnerVectorSizes_ = PHX::View<size_t*>(
"tangentInnerVectorSizes_",gatherFields_.size());
290 auto tangentInnerVectorSizes_host = Kokkos::create_mirror_view(tangentInnerVectorSizes_);
291 size_t inner_vector_max_size = 0;
292 for (std::size_t fd = 0; fd < tangentFields_.size(); ++fd) {
293 inner_vector_max_size = std::max(inner_vector_max_size,tangentFields_[fd].size());
294 tangentInnerVectorSizes_host(fd) = tangentFields_[fd].size();
296 Kokkos::deep_copy(tangentInnerVectorSizes_,tangentInnerVectorSizes_host);
298 gatherFieldsVoV_.initialize(
"GatherSolution_Teptra<Tangent>::gatherFieldsVoV_",gatherFields_.size());
299 tangentFieldsVoV_.initialize(
"GatherSolution_Teptra<Tangent>::tangentFieldsVoV_",gatherFields_.size(),inner_vector_max_size);
301 for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
302 const std::string& fieldName = indexerNames_[fd];
303 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
304 gatherFieldsVoV_.addView(gatherFields_[fd].get_static_view(),fd);
306 if (has_tangent_fields_) {
307 for (std::size_t i=0; i<tangentFields_[fd].size(); ++i) {
308 tangentFieldsVoV_.addView(tangentFields_[fd][i].get_static_view(),fd,i);
313 gatherFieldsVoV_.syncHostToDevice();
314 tangentFieldsVoV_.syncHostToDevice();
316 indexerNames_.clear();
320 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
327 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(d.gedc->getDataObject(globalDataKey_));
329 if(tpetraContainer_==Teuchos::null) {
332 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(loc);
337 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
344 std::string blockId = this->wda(workset).block_id;
347 if (useTimeDerivativeSolutionVector_)
348 x = tpetraContainer_->get_dxdt();
350 x = tpetraContainer_->get_x();
353 auto cellLocalIdsKokkos = this->wda(workset).getLocalCellIDs();
354 auto lids = globalIndexer_->getLIDs();
355 auto gidFieldOffsetsVoV = Teuchos::rcp_dynamic_cast<
const panzer::DOFManager>(globalIndexer_,
true)->getGIDFieldOffsetsKokkos(blockId,fieldIds_);
356 auto gidFieldOffsets = gidFieldOffsetsVoV.getViewDevice();
357 auto gatherFieldsDevice = gatherFieldsVoV_.getViewDevice();
358 auto x_view = x->getLocalViewDevice(Tpetra::Access::ReadOnly);
359 auto tangentInnerVectorSizes = this->tangentInnerVectorSizes_;
361 if (has_tangent_fields_) {
362 auto tangentFieldsDevice = tangentFieldsVoV_.getViewDevice();
363 Kokkos::parallel_for(
"GatherSolutionTpetra<Tangent>",cellLocalIdsKokkos.extent(0),KOKKOS_LAMBDA(
const int worksetCellIndex) {
364 for (
size_t fieldIndex = 0; fieldIndex < gidFieldOffsets.extent(0); ++fieldIndex) {
365 for(
size_t basis=0;basis<gidFieldOffsets(fieldIndex).extent(0);basis++) {
366 int offset = gidFieldOffsets(fieldIndex)(basis);
367 LO lid =
lids(cellLocalIdsKokkos(worksetCellIndex),offset);
368 auto gf_ref = (gatherFieldsDevice[fieldIndex])(worksetCellIndex,basis);
369 gf_ref.val() = x_view(lid,0);
370 for (std::size_t i=0; i<tangentInnerVectorSizes(fieldIndex); ++i) {
371 gf_ref.fastAccessDx(i) = tangentFieldsDevice(fieldIndex,i)(worksetCellIndex,basis);
378 Kokkos::parallel_for(
"GatherSolutionTpetra<Tangent>",cellLocalIdsKokkos.extent(0),KOKKOS_LAMBDA(
const int worksetCellIndex) {
379 for (
size_t fieldIndex = 0; fieldIndex < gidFieldOffsets.extent(0); ++fieldIndex) {
380 for(
size_t basis=0;basis<gidFieldOffsets(fieldIndex).extent(0);basis++) {
381 int offset = gidFieldOffsets(fieldIndex)(basis);
382 LO lid =
lids(cellLocalIdsKokkos(worksetCellIndex),offset);
383 reference_type gf_ref = (gatherFieldsDevice[fieldIndex])(worksetCellIndex,basis);
384 gf_ref.val() = x_view(lid,0);
395 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
400 : globalIndexer_(indexer)
407 const std::vector<std::string> & names = input.
getDofNames();
419 gatherFields_.resize(names.size());
420 scratch_offsets_.resize(names.size());
421 for (std::size_t fd = 0; fd < names.size(); ++fd) {
423 gatherFields_[fd] = f;
424 this->addEvaluatedField(gatherFields_[fd]);
427 this->addUnsharedField(gatherFields_[fd].fieldTag().clone());
431 std::string firstName =
"<none>";
433 firstName = names[0];
436 if(disableSensitivities_) {
437 std::string n =
"GatherSolution (Tpetra, No Sensitivities): "+firstName+
" (Jacobian)";
441 std::string n =
"GatherSolution (Tpetra): "+firstName+
" (Jacobian) ";
447 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
454 fieldIds_.resize(gatherFields_.size());
456 const Workset & workset_0 = (*d.worksets_)[0];
457 std::string blockId = this->wda(workset_0).
block_id;
459 for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
461 const std::string& fieldName = indexerNames_[fd];
462 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
464 int fieldNum = fieldIds_[fd];
465 const std::vector<int> &
offsets = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
466 scratch_offsets_[fd] = PHX::View<int*>(
"offsets",offsets.size());
470 scratch_lids_ = PHX::View<LO**>(
"lids",gatherFields_[0].extent(0),
471 globalIndexer_->getElementBlockGIDCount(blockId));
473 indexerNames_.clear();
477 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
483 using Teuchos::rcp_dynamic_cast;
490 if(!disableSensitivities_) {
491 if(d.first_sensitivities_name==sensitivitiesName_)
492 applySensitivities_ =
true;
494 applySensitivities_ =
false;
497 applySensitivities_ =
false;
504 std::string post = useTimeDerivativeSolutionVector_ ?
" - Xdot" :
" - X";
505 if(d.gedc->containsDataObject(globalDataKey_+post)) {
506 ged = d.gedc->getDataObject(globalDataKey_+post);
508 RCP<RO_GED> ro_ged = rcp_dynamic_cast<RO_GED>(ged,
true);
510 x_vector = ro_ged->getGhostedVector_Tpetra();
515 ged = d.gedc->getDataObject(globalDataKey_);
519 RCP<LOC> tpetraContainer = rcp_dynamic_cast<LOC>(ged);
522 if(loc_pair!=Teuchos::null) {
525 tpetraContainer = rcp_dynamic_cast<LOC>(loc);
528 if(tpetraContainer!=Teuchos::null) {
529 if (useTimeDerivativeSolutionVector_)
530 x_vector = tpetraContainer->get_dxdt();
532 x_vector = tpetraContainer->get_x();
540 RCP<RO_GED> ro_ged = rcp_dynamic_cast<RO_GED>(ged,
true);
542 x_vector = ro_ged->getGhostedVector_Tpetra();
547 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
552 std::string blockId = this->wda(workset).block_id;
554 double seed_value = 0.0;
555 if (useTimeDerivativeSolutionVector_) {
556 seed_value = workset.alpha;
558 else if (gatherSeedIndex_<0) {
559 seed_value = workset.beta;
561 else if(!useTimeDerivativeSolutionVector_) {
562 seed_value = workset.gather_seeds[gatherSeedIndex_];
571 if(!applySensitivities_)
577 functor_data.dos = 0;
578 if (this->wda.getDetailsIndex() == 1)
581 functor_data.dos = globalIndexer_->getElementBlockGIDCount(workset.details(0).block_id);
585 bool use_seed =
true;
589 globalIndexer_->getElementLIDs(this->wda(workset).cell_local_ids_k,scratch_lids_);
594 functor_data.x_data = x_vector->getLocalViewDevice(Tpetra::Access::ReadOnly);
595 functor_data.seed_value = seed_value;
596 functor_data.lids = scratch_lids_;
599 for(std::size_t fieldIndex=0;
600 fieldIndex<gatherFields_.size();fieldIndex++) {
603 functor_data.offsets = scratch_offsets_[fieldIndex];
604 functor_data.field = gatherFields_[fieldIndex];
607 Kokkos::parallel_for(workset.num_cells,*
this);
609 Kokkos::parallel_for(Kokkos::RangePolicy<PHX::Device,NoSeed>(0,workset.num_cells),*
this);
615 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
616 KOKKOS_INLINE_FUNCTION
621 for(std::size_t basis=0;basis<functor_data.offsets.extent(0);basis++) {
622 int offset = functor_data.offsets(basis);
623 LO lid = functor_data.lids(worksetCellIndex,offset);
626 if (functor_data.dos == 0)
627 functor_data.field(worksetCellIndex,basis).val() = functor_data.x_data(lid,0);
629 functor_data.field(worksetCellIndex,basis) =
ScalarT(functor_data.x_data(lid,0));
631 functor_data.field(worksetCellIndex,basis).fastAccessDx(functor_data.dos + offset) = functor_data.seed_value;
636 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
637 KOKKOS_INLINE_FUNCTION
642 for(std::size_t basis=0;basis<functor_data.offsets.extent(0);basis++) {
643 int offset = functor_data.offsets(basis);
644 LO lid = functor_data.lids(worksetCellIndex,offset);
647 functor_data.field(worksetCellIndex,basis).val() = functor_data.x_data(lid,0);
panzer::Traits::Jacobian::ScalarT ScalarT
PHX::View< const int * > offsets
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...
std::string block_id
DEPRECATED - use: getElementBlock()
#define TEUCHOS_ASSERT(assertion_test)
Teuchos::RCP< PHX::DataLayout > functional
<Cell,Basis> or <Cell,Basis>
Kokkos::View< const LO **, Kokkos::LayoutRight, PHX::Device > lids