Panzer  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
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<LO,GO> > & 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 
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<LO,GO> > & 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();
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] = Kokkos::View<int*,PHX::Device>("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  PHX::Device::fence();
173 
174  maxElementBlockGIDCount = std::max(fieldGlobalIndexers_[fd]->getElementBlockGIDCount(blockId),maxElementBlockGIDCount);
175  }
176 
177  // We will use one workset lid view for all fields, but has to be
178  // sized big enough to hold the largest elementBlockGIDCount in the
179  // ProductVector.
180  worksetLIDs_ = Kokkos::View<LO**,PHX::Device>("GatherSolution_BlockedTpetra(Residual):worksetLIDs",
181  gatherFields_[0].extent(0),
182  maxElementBlockGIDCount);
183 
184  indexerNames_.clear(); // Don't need this anymore
185 }
186 
187 // **********************************************************************
188 template <typename TRAITS,typename S,typename LO,typename GO,typename NodeT>
190 preEvaluate(typename TRAITS::PreEvalData d)
191 {
192  // extract linear object container
193  blockedContainer_ = Teuchos::rcp_dynamic_cast<const ContainerType>(d.gedc->getDataObject(globalDataKey_),true);
194 }
195 
196 // **********************************************************************
197 template <typename TRAITS,typename S,typename LO,typename GO,typename NodeT>
199 evaluateFields(typename TRAITS::EvalData workset)
200 {
201  using Teuchos::RCP;
202  using Teuchos::rcp_dynamic_cast;
203  using Thyra::VectorBase;
205 
206  const auto& localCellIds = this->wda(workset).cell_local_ids_k;
207 
208  RCP<ProductVectorBase<ScalarT>> thyraBlockSolution;
209  if (useTimeDerivativeSolutionVector_)
210  thyraBlockSolution = rcp_dynamic_cast<ProductVectorBase<ScalarT>>(blockedContainer_->get_dxdt(),true);
211  else
212  thyraBlockSolution = rcp_dynamic_cast<ProductVectorBase<ScalarT>>(blockedContainer_->get_x(),true);
213 
214  // Loop over gathered fields
215  int currentWorksetLIDSubBlock = -1;
216  for (std::size_t fieldIndex = 0; fieldIndex < gatherFields_.size(); fieldIndex++) {
217  // workset LIDs only change for different sub blocks
218  if (productVectorBlockIndex_[fieldIndex] != currentWorksetLIDSubBlock) {
219  fieldGlobalIndexers_[fieldIndex]->getElementLIDs(localCellIds,worksetLIDs_);
220  currentWorksetLIDSubBlock = productVectorBlockIndex_[fieldIndex];
221  }
222 
223  const auto& tpetraSolution = *((rcp_dynamic_cast<Thyra::TpetraVector<ScalarT,LO,GO,NodeT>>(thyraBlockSolution->getNonconstVectorBlock(productVectorBlockIndex_[fieldIndex]),true))->getTpetraVector());
224  const auto& kokkosSolution = tpetraSolution.template getLocalView<PHX::mem_space>();
225 
226  // Class data fields for lambda capture
227  const auto& fieldOffsets = fieldOffsets_[fieldIndex];
228  const auto& worksetLIDs = worksetLIDs_;
229  const auto& fieldValues = gatherFields_[fieldIndex];
230 
231  Kokkos::parallel_for(Kokkos::RangePolicy<PHX::Device>(0,workset.num_cells), KOKKOS_LAMBDA (const int& cell) {
232  for(int basis=0; basis < static_cast<int>(fieldOffsets.size()); ++basis) {
233  const int lid = worksetLIDs(cell,fieldOffsets(basis));
234  fieldValues(cell,basis) = kokkosSolution(lid,0);
235  }
236  });
237  }
238 
239 }
240 
241 // **********************************************************************
242 // Specialization: Tangent
243 // **********************************************************************
244 
245 template <typename TRAITS,typename S,typename LO,typename GO,typename NodeT>
248  const Teuchos::RCP<const BlockedDOFManager<LO,GO> > & indexer,
249  const Teuchos::ParameterList& p)
250  : gidIndexer_(indexer)
251  , has_tangent_fields_(false)
252 {
253  typedef std::vector< std::vector<std::string> > vvstring;
254 
255  GatherSolution_Input input;
256  input.setParameterList(p);
257 
258  const std::vector<std::string> & names = input.getDofNames();
260  const vvstring & tangent_field_names = input.getTangentNames();
261 
262  indexerNames_ = input.getIndexerNames();
263  useTimeDerivativeSolutionVector_ = input.useTimeDerivativeSolutionVector();
264  globalDataKey_ = input.getGlobalDataKey();
265 
266  // allocate fields
267  gatherFields_.resize(names.size());
268  for (std::size_t fd = 0; fd < names.size(); ++fd) {
269  gatherFields_[fd] =
270  PHX::MDField<ScalarT,Cell,NODE>(names[fd],basis->functional);
271  this->addEvaluatedField(gatherFields_[fd]);
272  }
273 
274  // Setup dependent tangent fields if requested
275  if (tangent_field_names.size()>0) {
276  TEUCHOS_ASSERT(gatherFields_.size() == tangent_field_names.size());
277 
278  has_tangent_fields_ = true;
279  tangentFields_.resize(gatherFields_.size());
280  for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
281  tangentFields_[fd].resize(tangent_field_names[fd].size());
282  for (std::size_t i=0; i<tangent_field_names[fd].size(); ++i) {
283  tangentFields_[fd][i] =
284  PHX::MDField<const ScalarT,Cell,NODE>(tangent_field_names[fd][i],basis->functional);
285  this->addDependentField(tangentFields_[fd][i]);
286  }
287  }
288  }
289 
290  // figure out what the first active name is
291  std::string firstName = "<none>";
292  if(names.size()>0)
293  firstName = names[0];
294 
295  std::string n = "GatherSolution (BlockedTpetra): "+firstName+" (Tangent)";
296  this->setName(n);
297 }
298 
299 // **********************************************************************
300 template <typename TRAITS,typename S,typename LO,typename GO,typename NodeT>
302 postRegistrationSetup(typename TRAITS::SetupData /* d */,
303  PHX::FieldManager<TRAITS>& /* fm */)
304 {
305  TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size());
306 
307  fieldIds_.resize(gatherFields_.size());
308 
309  for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
310  // get field ID from DOF manager
311  const std::string& fieldName = indexerNames_[fd];
312  fieldIds_[fd] = gidIndexer_->getFieldNum(fieldName);
313  }
314 
315  indexerNames_.clear(); // Don't need this anymore
316 }
317 
318 // **********************************************************************
319 template <typename TRAITS,typename S,typename LO,typename GO,typename NodeT>
321 preEvaluate(typename TRAITS::PreEvalData d)
322 {
323  // extract linear object container
324  blockedContainer_ = Teuchos::rcp_dynamic_cast<const ContainerType>(d.gedc->getDataObject(globalDataKey_),true);
325 }
326 
327 // **********************************************************************
328 template <typename TRAITS,typename S,typename LO,typename GO,typename NodeT>
330 evaluateFields(typename TRAITS::EvalData workset)
331 {
332  using Teuchos::RCP;
333  using Teuchos::ArrayRCP;
334  using Teuchos::ptrFromRef;
335  using Teuchos::rcp_dynamic_cast;
336 
337  using Thyra::VectorBase;
338  using Thyra::SpmdVectorBase;
340 
341  Teuchos::FancyOStream out(Teuchos::rcpFromRef(std::cout));
342  out.setShowProcRank(true);
343  out.setOutputToRootOnly(-1);
344 
345  std::vector<std::pair<int,GO> > GIDs;
346  std::vector<LO> LIDs;
347 
348  // for convenience pull out some objects from workset
349  std::string blockId = this->wda(workset).block_id;
350  const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
351 
353  if (useTimeDerivativeSolutionVector_)
354  x = rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer_->get_dxdt());
355  else
356  x = rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer_->get_x());
357 
358  // gather operation for each cell in workset
359  for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
360  LO cellLocalId = localCellIds[worksetCellIndex];
361 
362  gidIndexer_->getElementGIDs(cellLocalId,GIDs,blockId);
363 
364  // caculate the local IDs for this element
365  LIDs.resize(GIDs.size());
366  for(std::size_t i=0;i<GIDs.size();i++) {
367  // used for doing local ID lookups
368  RCP<const MapType> x_map = blockedContainer_->getMapForBlock(GIDs[i].first);
369 
370  LIDs[i] = x_map->getLocalElement(GIDs[i].second);
371  }
372 
373  // loop over the fields to be gathered
375  for (std::size_t fieldIndex=0; fieldIndex<gatherFields_.size();fieldIndex++) {
376  int fieldNum = fieldIds_[fieldIndex];
377  int indexerId = gidIndexer_->getFieldBlock(fieldNum);
378 
379  // grab local data for inputing
380  RCP<SpmdVectorBase<double> > block_x = rcp_dynamic_cast<SpmdVectorBase<double> >(x->getNonconstVectorBlock(indexerId));
381  block_x->getLocalData(ptrFromRef(local_x));
382 
383  const std::vector<int> & elmtOffset = gidIndexer_->getGIDFieldOffsets(blockId,fieldNum);
384 
385  // loop over basis functions and fill the fields
386  for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
387  int offset = elmtOffset[basis];
388  int lid = LIDs[offset];
389 
390  if (!has_tangent_fields_)
391  (gatherFields_[fieldIndex])(worksetCellIndex,basis) = local_x[lid];
392  else {
393  (gatherFields_[fieldIndex])(worksetCellIndex,basis).val() = local_x[lid];
394  for (std::size_t i=0; i<tangentFields_[fieldIndex].size(); ++i)
395  (gatherFields_[fieldIndex])(worksetCellIndex,basis).fastAccessDx(i) =
396  tangentFields_[fieldIndex][i](worksetCellIndex,basis).val();
397  }
398  }
399  }
400  }
401 }
402 
403 // **********************************************************************
404 // Specialization: Jacobian
405 // **********************************************************************
406 
407 template <typename TRAITS,typename S,typename LO,typename GO,typename NodeT>
410  const Teuchos::RCP<const BlockedDOFManager<LO,GO> > & indexer,
411  const Teuchos::ParameterList& p)
412  : globalIndexer_(indexer)
413 {
414  GatherSolution_Input input;
415  input.setParameterList(p);
416 
417  const std::vector<std::string> & names = input.getDofNames();
419 
420  indexerNames_ = input.getIndexerNames();
421  useTimeDerivativeSolutionVector_ = input.useTimeDerivativeSolutionVector();
422  globalDataKey_ = input.getGlobalDataKey();
423 
424  disableSensitivities_ = !input.firstSensitivitiesAvailable();
425 
426  gatherFields_.resize(names.size());
427  for (std::size_t fd = 0; fd < names.size(); ++fd) {
428  PHX::MDField<ScalarT,Cell,NODE> f(names[fd],basis->functional);
429  gatherFields_[fd] = f;
430  this->addEvaluatedField(gatherFields_[fd]);
431  }
432 
433  // figure out what the first active name is
434  std::string firstName = "<none>";
435  if(names.size()>0)
436  firstName = names[0];
437 
438  // print out convenience
439  if(disableSensitivities_) {
440  std::string n = "GatherSolution (BlockedTpetra, No Sensitivities): "+firstName+" (Jacobian)";
441  this->setName(n);
442  }
443  else {
444  std::string n = "GatherSolution (BlockedTpetra): "+firstName+" ("+PHX::typeAsString<EvalT>()+") ";
445  this->setName(n);
446  }
447 }
448 
449 // **********************************************************************
450 template <typename TRAITS,typename S,typename LO,typename GO,typename NodeT>
452 postRegistrationSetup(typename TRAITS::SetupData d,
453  PHX::FieldManager<TRAITS>& /* fm */)
454 {
455  TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size());
456 
457  const Workset & workset_0 = (*d.worksets_)[0];
458  const std::string blockId = this->wda(workset_0).block_id;
459 
460  fieldIds_.resize(gatherFields_.size());
461  fieldOffsets_.resize(gatherFields_.size());
462  productVectorBlockIndex_.resize(gatherFields_.size());
463  int maxElementBlockGIDCount = -1;
464  for (std::size_t fd = 0; fd < gatherFields_.size(); ++fd) {
465 
466  const std::string fieldName = indexerNames_[fd];
467  const int globalFieldNum = globalIndexer_->getFieldNum(fieldName); // Field number in the aggregate BlockDOFManager
468  productVectorBlockIndex_[fd] = globalIndexer_->getFieldBlock(globalFieldNum);
469  const auto& subGlobalIndexer = globalIndexer_->getFieldDOFManagers()[productVectorBlockIndex_[fd]];
470  fieldIds_[fd] = subGlobalIndexer->getFieldNum(fieldName); // Field number in the sub-global-indexer
471 
472  const std::vector<int>& offsets = subGlobalIndexer->getGIDFieldOffsets(blockId,fieldIds_[fd]);
473  fieldOffsets_[fd] = Kokkos::View<int*,PHX::Device>("GatherSolution_BlockedTpetra(Jacobian):fieldOffsets",offsets.size());
474  auto hostOffsets = Kokkos::create_mirror_view(fieldOffsets_[fd]);
475  for (std::size_t i=0; i < offsets.size(); ++i)
476  hostOffsets(i) = offsets[i];
477  Kokkos::deep_copy(fieldOffsets_[fd], hostOffsets);
478  maxElementBlockGIDCount = std::max(subGlobalIndexer->getElementBlockGIDCount(blockId),maxElementBlockGIDCount);
479  PHX::Device::fence();
480  }
481 
482  // We will use one workset lid view for all fields, but has to be
483  // sized big enough to hold the largest elementBlockGIDCount in the
484  // ProductVector.
485  worksetLIDs_ = Kokkos::View<LO**,PHX::Device>("ScatterResidual_BlockedTpetra(Residual):worksetLIDs",
486  gatherFields_[0].extent(0),
487  maxElementBlockGIDCount);
488 
489  // Compute the block offsets
490  const auto& blockGlobalIndexers = globalIndexer_->getFieldDOFManagers();
491  const int numBlocks = static_cast<int>(globalIndexer_->getFieldDOFManagers().size());
492  blockOffsets_ = Kokkos::View<LO*,PHX::Device>("GatherSolution_BlockedTpetra(Jacobian):blockOffsets_",
493  numBlocks+1); // Number of blocks, plus a sentinel
494  const auto hostBlockOffsets = Kokkos::create_mirror_view(blockOffsets_);
495  for (int blk=0;blk<numBlocks;++blk) {
496  int blockOffset = globalIndexer_->getBlockGIDOffset(blockId,blk);
497  hostBlockOffsets(blk) = blockOffset;
498  }
499  blockOffsets_(numBlocks) = blockOffsets_(numBlocks-1) + blockGlobalIndexers[blockGlobalIndexers.size()-1]->getElementBlockGIDCount(blockId);
500  Kokkos::deep_copy(blockOffsets_,hostBlockOffsets);
501 
502  indexerNames_.clear(); // Don't need this anymore
503  PHX::Device::fence();
504 }
505 
506 template <typename TRAITS,typename S,typename LO,typename GO,typename NodeT>
508 preEvaluate(typename TRAITS::PreEvalData d)
509 {
510  // extract linear object container
511  blockedContainer_ = Teuchos::rcp_dynamic_cast<const ContainerType>(d.gedc->getDataObject(globalDataKey_),true);
512 }
513 
514 // **********************************************************************
515 template <typename TRAITS,typename S,typename LO,typename GO,typename NodeT>
517 evaluateFields(typename TRAITS::EvalData workset)
518 {
519  using Teuchos::RCP;
520  using Teuchos::rcp_dynamic_cast;
521  using Thyra::VectorBase;
523 
524  const auto& localCellIds = this->wda(workset).cell_local_ids_k;
525 
526  RealType seedValue = RealType(0.0);
527  RCP<ProductVectorBase<double>> blockedSolution;
528  if (useTimeDerivativeSolutionVector_) {
529  blockedSolution = rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer_->get_dxdt());
530  seedValue = workset.alpha;
531  }
532  else {
533  blockedSolution = rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer_->get_x());
534  seedValue = workset.beta;
535  }
536 
537  // turn off sensitivies: this may be faster if we don't expand the term
538  // but I suspect not because anywhere it is used the full complement of
539  // sensitivies will be needed anyway.
540  if(disableSensitivities_)
541  seedValue = 0.0;
542 
543  const int numFieldBlocks = globalIndexer_->getNumFieldBlocks();
544 
545  // Loop over fields to gather
546  int currentWorksetLIDSubBlock = -1;
547  for (std::size_t fieldIndex = 0; fieldIndex < gatherFields_.size(); fieldIndex++) {
548  // workset LIDs only change if in different sub blocks
549  if (productVectorBlockIndex_[fieldIndex] != currentWorksetLIDSubBlock) {
550  const auto& blockIndexer = globalIndexer_->getFieldDOFManagers()[productVectorBlockIndex_[fieldIndex]];
551  blockIndexer->getElementLIDs(localCellIds,worksetLIDs_);
552  currentWorksetLIDSubBlock = productVectorBlockIndex_[fieldIndex];
553  }
554 
555  const int blockRowIndex = productVectorBlockIndex_[fieldIndex];
556  const auto& subblockSolution = *((rcp_dynamic_cast<Thyra::TpetraVector<RealType,LO,GO,NodeT>>(blockedSolution->getNonconstVectorBlock(blockRowIndex),true))->getTpetraVector());
557  const auto kokkosSolution = subblockSolution.template getLocalView<PHX::mem_space>();
558 
559  // Class data fields for lambda capture
560  const Kokkos::View<const int*,PHX::Device> fieldOffsets = fieldOffsets_[fieldIndex];
561  const Kokkos::View<const LO**,PHX::Device> worksetLIDs = worksetLIDs_;
562  const PHX::View<ScalarT**> fieldValues = gatherFields_[fieldIndex].get_static_view();
563  const Kokkos::View<const LO*,PHX::Device> blockOffsets = blockOffsets_;
564  const int blockStart = blockOffsets(blockRowIndex);
565  const int numDerivatives = blockOffsets(numFieldBlocks);
566 
567  Kokkos::parallel_for(Kokkos::RangePolicy<PHX::Device>(0,workset.num_cells), KOKKOS_LAMBDA (const int& cell) {
568  for (int basis=0; basis < static_cast<int>(fieldOffsets.size()); ++basis) {
569  const int rowLID = worksetLIDs(cell,fieldOffsets(basis));
570  fieldValues(cell,basis) = ScalarT(numDerivatives,kokkosSolution(rowLID,0));
571  fieldValues(cell,basis).fastAccessDx(blockStart+fieldOffsets(basis)) = seedValue;
572  }
573  });
574 
575  }
576 }
577 
578 // **********************************************************************
579 
580 #endif
const std::vector< std::vector< std::string > > & getTangentNames() const
Get the name of the tangent fields (tangent only)
const std::vector< std::string > & getDofNames() const
The names of the DOFs to be gathered (all types)
void setParameterList(const Teuchos::ParameterList &pl)
GatherSolution_BlockedTpetra(const Teuchos::RCP< const BlockedDOFManager< LO, GO > > &indexer)
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &vm)
basic_FancyOStream & setShowProcRank(const bool showProcRank)
T & get(const std::string &name, T def_value)
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)
Gathers solution values from the Newton solution vector into the nodal fields of the field manager...
bool firstSensitivitiesAvailable()
Are first derivative sensitivities enabled or disabled? (Jacobian and Hessian)
bool useTimeDerivativeSolutionVector() const
Gather a time derivative vector? (all types)
basic_FancyOStream & setOutputToRootOnly(const int rootRank)
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
#define TEUCHOS_ASSERT(assertion_test)
Teuchos::RCP< PHX::DataLayout > functional
&lt;Cell,Basis&gt; or &lt;Cell,Basis&gt;
Kokkos::View< const int *, PHX::Device > offsets