Panzer  Version of the Day
Panzer_GatherSolution_BlockedTpetra_impl.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Panzer: A partial differential equation assembly
5 // engine for strongly coupled complex multiphysics systems
6 // Copyright (2011) Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Roger P. Pawlowski (rppawlo@sandia.gov) and
39 // Eric C. Cyr (eccyr@sandia.gov)
40 // ***********************************************************************
41 // @HEADER
42 
43 #ifndef PANZER_GATHER_SOLUTION_BLOCKED_EPETRA_IMPL_HPP
44 #define PANZER_GATHER_SOLUTION_BLOCKED_EPETRA_IMPL_HPP
45 
46 #include "Teuchos_Assert.hpp"
47 #include "Phalanx_DataLayout.hpp"
48 
49 #include "Panzer_GlobalIndexer.hpp"
51 #include "Panzer_PureBasis.hpp"
52 #include "Panzer_TpetraLinearObjFactory.hpp"
56 
57 #include "Teuchos_FancyOStream.hpp"
58 
59 #include "Thyra_SpmdVectorBase.hpp"
60 #include "Thyra_ProductVectorBase.hpp"
61 
62 #include "Tpetra_Map.hpp"
63 
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)
69 {
70  const std::vector<std::string>& names =
71  *(p.get< Teuchos::RCP< std::vector<std::string> > >("DOF Names"));
72 
73  Teuchos::RCP<panzer::PureBasis> basis =
74  p.get< Teuchos::RCP<panzer::PureBasis> >("Basis");
75 
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());
79  }
80 
81  this->setName("Gather Solution");
82 }
83 
84 // **********************************************************************
85 // Specialization: Residual
86 // **********************************************************************
87 
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)
95 {
96  typedef std::vector< std::vector<std::string> > vvstring;
97 
99  input.setParameterList(p);
100 
101  const std::vector<std::string> & names = input.getDofNames();
102  Teuchos::RCP<const panzer::PureBasis> basis = input.getBasis();
103  const vvstring & tangent_field_names = input.getTangentNames();
104 
105  indexerNames_ = input.getIndexerNames();
106  useTimeDerivativeSolutionVector_ = input.useTimeDerivativeSolutionVector();
107  globalDataKey_ = input.getGlobalDataKey();
108 
109  // allocate fields
110  gatherFields_.resize(names.size());
111  for (std::size_t fd = 0; fd < names.size(); ++fd) {
112  gatherFields_[fd] =
113  PHX::MDField<ScalarT,Cell,NODE>(names[fd],basis->functional);
114  this->addEvaluatedField(gatherFields_[fd]);
115  }
116 
117  // Setup dependent tangent fields if requested
118  if (tangent_field_names.size()>0) {
119  TEUCHOS_ASSERT(gatherFields_.size() == tangent_field_names.size());
120 
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]);
129  }
130  }
131  }
132 
133  // figure out what the first active name is
134  std::string firstName = "<none>";
135  if(names.size()>0)
136  firstName = names[0];
137 
138  std::string n = "GatherSolution (BlockedTpetra): "+firstName+" (Residual)";
139  this->setName(n);
140 }
141 
142 // **********************************************************************
143 template <typename TRAITS,typename S,typename LO,typename GO,typename NodeT>
145 postRegistrationSetup(typename TRAITS::SetupData d,
146  PHX::FieldManager<TRAITS>& /* fm */)
147 {
148  TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size());
149 
150  const Workset & workset_0 = (*d.worksets_)[0];
151  const std::string blockId = this->wda(workset_0).block_id;
152 
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) {
159  // get field ID from DOF manager
160  const std::string& fieldName = indexerNames_[fd];
161  const int globalFieldNum = globalIndexer_->getFieldNum(fieldName); // Field number in the aggregate BlockDOFManager
162  productVectorBlockIndex_[fd] = globalIndexer_->getFieldBlock(globalFieldNum);
163  fieldGlobalIndexers_[fd] = globalIndexer_->getFieldDOFManagers()[productVectorBlockIndex_[fd]];
164  fieldIds_[fd] = fieldGlobalIndexers_[fd]->getFieldNum(fieldName); // Field number in the sub-global-indexer
165 
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);
172 
173  maxElementBlockGIDCount = std::max(fieldGlobalIndexers_[fd]->getElementBlockGIDCount(blockId),maxElementBlockGIDCount);
174  }
175 
176  // We will use one workset lid view for all fields, but has to be
177  // sized big enough to hold the largest elementBlockGIDCount in the
178  // ProductVector.
179  worksetLIDs_ = PHX::View<LO**>("GatherSolution_BlockedTpetra(Residual):worksetLIDs",
180  gatherFields_[0].extent(0),
181  maxElementBlockGIDCount);
182 
183  indexerNames_.clear(); // Don't need this anymore
184 }
185 
186 // **********************************************************************
187 template <typename TRAITS,typename S,typename LO,typename GO,typename NodeT>
189 preEvaluate(typename TRAITS::PreEvalData d)
190 {
191  // extract linear object container
192  blockedContainer_ = Teuchos::rcp_dynamic_cast<const ContainerType>(d.gedc->getDataObject(globalDataKey_),true);
193 }
194 
195 // **********************************************************************
196 template <typename TRAITS,typename S,typename LO,typename GO,typename NodeT>
198 evaluateFields(typename TRAITS::EvalData workset)
199 {
200  using Teuchos::RCP;
201  using Teuchos::rcp_dynamic_cast;
202  using Thyra::VectorBase;
204 
205  const PHX::View<const int*>& localCellIds = this->wda(workset).cell_local_ids_k;
206 
207  RCP<ProductVectorBase<ScalarT>> thyraBlockSolution;
208  if (useTimeDerivativeSolutionVector_)
209  thyraBlockSolution = rcp_dynamic_cast<ProductVectorBase<ScalarT>>(blockedContainer_->get_dxdt(),true);
210  else
211  thyraBlockSolution = rcp_dynamic_cast<ProductVectorBase<ScalarT>>(blockedContainer_->get_x(),true);
212 
213  // Loop over gathered fields
214  int currentWorksetLIDSubBlock = -1;
215  for (std::size_t fieldIndex = 0; fieldIndex < gatherFields_.size(); fieldIndex++) {
216  // workset LIDs only change for different sub blocks
217  if (productVectorBlockIndex_[fieldIndex] != currentWorksetLIDSubBlock) {
218  fieldGlobalIndexers_[fieldIndex]->getElementLIDs(localCellIds,worksetLIDs_);
219  currentWorksetLIDSubBlock = productVectorBlockIndex_[fieldIndex];
220  }
221 
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);
224 
225  // Class data fields for lambda capture
226  const auto& fieldOffsets = fieldOffsets_[fieldIndex];
227  const auto& worksetLIDs = worksetLIDs_;
228  const auto& fieldValues = gatherFields_[fieldIndex];
229 
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);
234  }
235  });
236  }
237 
238 }
239 
240 // **********************************************************************
241 // Specialization: Tangent
242 // **********************************************************************
243 
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)
251 {
252  typedef std::vector< std::vector<std::string> > vvstring;
253 
254  GatherSolution_Input input;
255  input.setParameterList(p);
256 
257  const std::vector<std::string> & names = input.getDofNames();
258  Teuchos::RCP<const panzer::PureBasis> basis = input.getBasis();
259  const vvstring & tangent_field_names = input.getTangentNames();
260 
261  indexerNames_ = input.getIndexerNames();
262  useTimeDerivativeSolutionVector_ = input.useTimeDerivativeSolutionVector();
263  globalDataKey_ = input.getGlobalDataKey();
264 
265  // allocate fields
266  gatherFields_.resize(names.size());
267  for (std::size_t fd = 0; fd < names.size(); ++fd) {
268  gatherFields_[fd] =
269  PHX::MDField<ScalarT,Cell,NODE>(names[fd],basis->functional);
270  this->addEvaluatedField(gatherFields_[fd]);
271  }
272 
273  // Setup dependent tangent fields if requested
274  if (tangent_field_names.size()>0) {
275  TEUCHOS_ASSERT(gatherFields_.size() == tangent_field_names.size());
276 
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]);
285  }
286  }
287  }
288 
289  // figure out what the first active name is
290  std::string firstName = "<none>";
291  if(names.size()>0)
292  firstName = names[0];
293 
294  std::string n = "GatherSolution (BlockedTpetra): "+firstName+" (Tangent)";
295  this->setName(n);
296 }
297 
298 // **********************************************************************
299 template <typename TRAITS,typename S,typename LO,typename GO,typename NodeT>
301 postRegistrationSetup(typename TRAITS::SetupData /* d */,
302  PHX::FieldManager<TRAITS>& /* fm */)
303 {
304  TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size());
305 
306  fieldIds_.resize(gatherFields_.size());
307 
308  for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
309  // get field ID from DOF manager
310  const std::string& fieldName = indexerNames_[fd];
311  fieldIds_[fd] = gidIndexer_->getFieldNum(fieldName);
312  }
313 
314  indexerNames_.clear(); // Don't need this anymore
315 }
316 
317 // **********************************************************************
318 template <typename TRAITS,typename S,typename LO,typename GO,typename NodeT>
320 preEvaluate(typename TRAITS::PreEvalData d)
321 {
322  // extract linear object container
323  blockedContainer_ = Teuchos::rcp_dynamic_cast<const ContainerType>(d.gedc->getDataObject(globalDataKey_),true);
324 }
325 
326 // **********************************************************************
327 template <typename TRAITS,typename S,typename LO,typename GO,typename NodeT>
329 evaluateFields(typename TRAITS::EvalData workset)
330 {
331  using Teuchos::RCP;
332  using Teuchos::ArrayRCP;
333  using Teuchos::ptrFromRef;
334  using Teuchos::rcp_dynamic_cast;
335 
336  using Thyra::VectorBase;
337  using Thyra::SpmdVectorBase;
339 
340  Teuchos::FancyOStream out(Teuchos::rcpFromRef(std::cout));
341  out.setShowProcRank(true);
342  out.setOutputToRootOnly(-1);
343 
344  std::vector<std::pair<int,GO> > GIDs;
345  std::vector<LO> LIDs;
346 
347  // for convenience pull out some objects from workset
348  std::string blockId = this->wda(workset).block_id;
349  const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
350 
351  Teuchos::RCP<ProductVectorBase<double> > x;
352  if (useTimeDerivativeSolutionVector_)
353  x = rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer_->get_dxdt());
354  else
355  x = rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer_->get_x());
356 
357  // gather operation for each cell in workset
358  for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
359  LO cellLocalId = localCellIds[worksetCellIndex];
360 
361  gidIndexer_->getElementGIDsPair(cellLocalId,GIDs,blockId);
362 
363  // caculate the local IDs for this element
364  LIDs.resize(GIDs.size());
365  for(std::size_t i=0;i<GIDs.size();i++) {
366  // used for doing local ID lookups
367  RCP<const MapType> x_map = blockedContainer_->getMapForBlock(GIDs[i].first);
368 
369  LIDs[i] = x_map->getLocalElement(GIDs[i].second);
370  }
371 
372  // loop over the fields to be gathered
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);
377 
378  // grab local data for inputing
379  RCP<SpmdVectorBase<double> > block_x = rcp_dynamic_cast<SpmdVectorBase<double> >(x->getNonconstVectorBlock(indexerId));
380  block_x->getLocalData(ptrFromRef(local_x));
381 
382  const std::vector<int> & elmtOffset = gidIndexer_->getGIDFieldOffsets(blockId,fieldNum);
383 
384  // loop over basis functions and fill the fields
385  for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
386  int offset = elmtOffset[basis];
387  int lid = LIDs[offset];
388 
389  if (!has_tangent_fields_)
390  (gatherFields_[fieldIndex])(worksetCellIndex,basis) = local_x[lid];
391  else {
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();
396  }
397  }
398  }
399  }
400 }
401 
402 // **********************************************************************
403 // Specialization: Jacobian
404 // **********************************************************************
405 
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)
412 {
413  GatherSolution_Input input;
414  input.setParameterList(p);
415 
416  const std::vector<std::string> & names = input.getDofNames();
417  Teuchos::RCP<const panzer::PureBasis> basis = input.getBasis();
418 
419  indexerNames_ = input.getIndexerNames();
420  useTimeDerivativeSolutionVector_ = input.useTimeDerivativeSolutionVector();
421  globalDataKey_ = input.getGlobalDataKey();
422 
423  disableSensitivities_ = !input.firstSensitivitiesAvailable();
424 
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]);
430  }
431 
432  // figure out what the first active name is
433  std::string firstName = "<none>";
434  if(names.size()>0)
435  firstName = names[0];
436 
437  // print out convenience
438  if(disableSensitivities_) {
439  std::string n = "GatherSolution (BlockedTpetra, No Sensitivities): "+firstName+" (Jacobian)";
440  this->setName(n);
441  }
442  else {
443  std::string n = "GatherSolution (BlockedTpetra): "+firstName+" ("+PHX::print<EvalT>()+") ";
444  this->setName(n);
445  }
446 }
447 
448 // **********************************************************************
449 template <typename TRAITS,typename S,typename LO,typename GO,typename NodeT>
451 postRegistrationSetup(typename TRAITS::SetupData d,
452  PHX::FieldManager<TRAITS>& /* fm */)
453 {
454  TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size());
455 
456  const Workset & workset_0 = (*d.worksets_)[0];
457  const std::string blockId = this->wda(workset_0).block_id;
458 
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) {
464 
465  const std::string fieldName = indexerNames_[fd];
466  const int globalFieldNum = globalIndexer_->getFieldNum(fieldName); // Field number in the aggregate BlockDOFManager
467  productVectorBlockIndex_[fd] = globalIndexer_->getFieldBlock(globalFieldNum);
468  const auto& subGlobalIndexer = globalIndexer_->getFieldDOFManagers()[productVectorBlockIndex_[fd]];
469  fieldIds_[fd] = subGlobalIndexer->getFieldNum(fieldName); // Field number in the sub-global-indexer
470 
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)
475  hostOffsets(i) = offsets[i];
476  Kokkos::deep_copy(fieldOffsets_[fd], hostOffsets);
477  maxElementBlockGIDCount = std::max(subGlobalIndexer->getElementBlockGIDCount(blockId),maxElementBlockGIDCount);
478  }
479 
480  // We will use one workset lid view for all fields, but has to be
481  // sized big enough to hold the largest elementBlockGIDCount in the
482  // ProductVector.
483  worksetLIDs_ = PHX::View<LO**>("ScatterResidual_BlockedTpetra(Residual):worksetLIDs",
484  gatherFields_[0].extent(0),
485  maxElementBlockGIDCount);
486 
487  // Compute the block offsets
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_",
491  numBlocks+1); // Number of blocks, plus a sentinel
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;
496  }
497  hostBlockOffsets(numBlocks) = hostBlockOffsets(numBlocks-1) + blockGlobalIndexers[blockGlobalIndexers.size()-1]->getElementBlockGIDCount(blockId);
498  Kokkos::deep_copy(blockOffsets_,hostBlockOffsets);
499 
500  indexerNames_.clear(); // Don't need this anymore
501 }
502 
503 template <typename TRAITS,typename S,typename LO,typename GO,typename NodeT>
505 preEvaluate(typename TRAITS::PreEvalData d)
506 {
507  // extract linear object container
508  blockedContainer_ = Teuchos::rcp_dynamic_cast<const ContainerType>(d.gedc->getDataObject(globalDataKey_),true);
509 }
510 
511 // **********************************************************************
512 template <typename TRAITS,typename S,typename LO,typename GO,typename NodeT>
514 evaluateFields(typename TRAITS::EvalData workset)
515 {
516  using Teuchos::RCP;
517  using Teuchos::rcp_dynamic_cast;
518  using Thyra::VectorBase;
520 
521  const auto& localCellIds = this->wda(workset).cell_local_ids_k;
522 
523  RealType seedValue = RealType(0.0);
524  RCP<ProductVectorBase<double>> blockedSolution;
525  if (useTimeDerivativeSolutionVector_) {
526  blockedSolution = rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer_->get_dxdt());
527  seedValue = workset.alpha;
528  }
529  else {
530  blockedSolution = rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer_->get_x());
531  seedValue = workset.beta;
532  }
533 
534  // turn off sensitivies: this may be faster if we don't expand the term
535  // but I suspect not because anywhere it is used the full complement of
536  // sensitivies will be needed anyway.
537  if(disableSensitivities_)
538  seedValue = 0.0;
539 
540  const int numFieldBlocks = globalIndexer_->getNumFieldBlocks();
541 
542  // Loop over fields to gather
543  int currentWorksetLIDSubBlock = -1;
544  for (std::size_t fieldIndex = 0; fieldIndex < gatherFields_.size(); fieldIndex++) {
545  // workset LIDs only change if in different sub blocks
546  if (productVectorBlockIndex_[fieldIndex] != currentWorksetLIDSubBlock) {
547  const auto& blockIndexer = globalIndexer_->getFieldDOFManagers()[productVectorBlockIndex_[fieldIndex]];
548  blockIndexer->getElementLIDs(localCellIds,worksetLIDs_);
549  currentWorksetLIDSubBlock = productVectorBlockIndex_[fieldIndex];
550  }
551 
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);
555 
556  // Class data fields for lambda capture
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);
565 
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;
571  }
572  });
573 
574  }
575 }
576 
577 // **********************************************************************
578 
579 #endif
std::string getGlobalDataKey() const
Name of the global evaluation data container to use for the source vector (all types) ...
Teuchos::RCP< const PureBasis > getBasis() const
Basis definiting layout of dof names (all types)
GatherSolution_BlockedTpetra(const Teuchos::RCP< const BlockedDOFManager > &indexer)
void setParameterList(const Teuchos::ParameterList &pl)
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &vm)
const std::vector< std::string > & getDofNames() const
The names of the DOFs to be gathered (all types)
PHX::View< const int * > offsets
Gathers solution values from the Newton solution vector into the nodal fields of the field manager...
bool useTimeDerivativeSolutionVector() const
Gather a time derivative vector? (all types)
bool firstSensitivitiesAvailable()
Are first derivative sensitivities enabled or disabled? (Jacobian and Hessian)
std::string block_id
DEPRECATED - use: getElementBlock()
PHX::MDField< ScalarT, panzer::Cell, panzer::BASIS > field
A field to which we&#39;ll contribute, or in which we&#39;ll store, the result of computing this integral...
const std::vector< std::string > & getIndexerNames() const
const std::vector< std::vector< std::string > > & getTangentNames() const
Get the name of the tangent fields (tangent only)