Panzer  Version of the Day
Panzer_GatherSolution_BlockedEpetra_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_GatherSolution_BlockedEpetra_impl_hpp__
44 #define __Panzer_GatherSolution_BlockedEpetra_impl_hpp__
45 
47 //
48 // Include Files
49 //
51 
52 // Panzer
56 #include "Panzer_PureBasis.hpp"
57 #include "Panzer_GlobalIndexer.hpp"
60 
61 // Phalanx
62 #include "Phalanx_DataLayout.hpp"
63 
64 // Teuchos
65 #include "Teuchos_Assert.hpp"
66 #include "Teuchos_FancyOStream.hpp"
67 
68 // Thyra
69 #include "Thyra_ProductVectorBase.hpp"
70 #include "Thyra_SpmdVectorBase.hpp"
71 
73 //
74 // Initializing Constructor (Residual Specialization)
75 //
77 template<typename TRAITS, typename LO, typename GO>
81  const std::vector<Teuchos::RCP<const GlobalIndexer>>&
82  indexers,
83  const Teuchos::ParameterList& p)
84  :
85  indexers_(indexers),
86  hasTangentFields_(false)
87 {
88  using panzer::PureBasis;
89  using PHX::MDField;
90  using PHX::print;
91  using std::size_t;
92  using std::string;
93  using std::vector;
94  using Teuchos::RCP;
95  using vvstring = std::vector<std::vector<std::string>>;
97  input.setParameterList(p);
98  const vector<string>& names = input.getDofNames();
99  RCP<const PureBasis> basis = input.getBasis();
100  const vvstring& tangentFieldNames = input.getTangentNames();
101  indexerNames_ = input.getIndexerNames();
102  useTimeDerivativeSolutionVector_ = input.useTimeDerivativeSolutionVector();
103  globalDataKey_ = input.getGlobalDataKey();
104 
105  // Allocate the fields.
106  int numFields(names.size());
107  gatherFields_.resize(numFields);
108  for (int fd(0); fd < numFields; ++fd)
109  {
110  gatherFields_[fd] =
111  MDField<ScalarT, Cell, NODE>(names[fd], basis->functional);
112  this->addEvaluatedField(gatherFields_[fd]);
113  } // end loop over names
114 
115  // Setup the dependent tangent fields, if requested.
116  if (tangentFieldNames.size() > 0)
117  {
118  TEUCHOS_ASSERT(gatherFields_.size() == tangentFieldNames.size());
119  hasTangentFields_ = true;
120  tangentFields_.resize(numFields);
121  for (int fd(0); fd < numFields; ++fd)
122  {
123  int numTangentFields(tangentFieldNames[fd].size());
124  tangentFields_[fd].resize(numTangentFields);
125  for (int i(0); i < numTangentFields; ++i)
126  {
127  tangentFields_[fd][i] =
128  MDField<const ScalarT, Cell, NODE>(tangentFieldNames[fd][i],
129  basis->functional);
130  this->addDependentField(tangentFields_[fd][i]);
131  } // end loop over tangentFieldNames[fd]
132  } // end loop over gatherFields_
133  } // end if we have tangent fields
134 
135  // Figure out what the first active name is.
136  string firstName("<none>");
137  if (numFields > 0)
138  firstName = names[0];
139  string n("GatherSolution (BlockedEpetra): " + firstName + " (" +
140  print<EvalT>() + ")");
141  this->setName(n);
142 } // end of Initializing Constructor (Residual Specialization)
143 
145 //
146 // postRegistrationSetup() (Residual Specialization)
147 //
149 template<typename TRAITS, typename LO, typename GO>
150 void
151 panzer::
154  typename TRAITS::SetupData /* d */,
155  PHX::FieldManager<TRAITS>& /* fm */)
156 {
157  using std::size_t;
158  using std::string;
159  TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size());
160  int numFields(gatherFields_.size());
161  indexerIds_.resize(numFields);
162  subFieldIds_.resize(numFields);
163  for (int fd(0); fd < numFields; ++fd)
164  {
165  // Get the field ID from the DOF manager.
166  const string& fieldName(indexerNames_[fd]);
167  indexerIds_[fd] = getFieldBlock(fieldName, indexers_);
168  subFieldIds_[fd] = indexers_[indexerIds_[fd]]->getFieldNum(fieldName);
169  TEUCHOS_ASSERT(indexerIds_[fd] >= 0);
170  } // end loop over gatherFields_
171  indexerNames_.clear();
172 } // end of postRegistrationSetup() (Residual Specialization)
173 
175 //
176 // preEvaluate() (Residual Specialization)
177 //
179 template<typename TRAITS, typename LO, typename GO>
180 void
181 panzer::
184  typename TRAITS::PreEvalData d)
185 {
186  using std::logic_error;
187  using std::string;
188  using Teuchos::RCP;
189  using Teuchos::rcp_dynamic_cast;
190  using Teuchos::typeName;
194  using GED = panzer::GlobalEvaluationData;
195 
196  // First try the refactored ReadOnly container.
197  RCP<GED> ged;
198  string post(useTimeDerivativeSolutionVector_ ? " - Xdot" : " - X");
199  if (d.gedc->containsDataObject(globalDataKey_ + post))
200  {
201  ged = d.gedc->getDataObject(globalDataKey_ + post);
202  xBvRoGed_ = rcp_dynamic_cast<BVROGED>(ged, true);
203  return;
204  } // end of the refactored ReadOnly way
205 
206  // Now try the old path.
207  {
208  ged = d.gedc->getDataObject(globalDataKey_);
209 
210  // Extract the linear object container.
211  auto roGed = rcp_dynamic_cast<const BVROGED>(ged);
212  auto beLoc = rcp_dynamic_cast<const BELOC>(ged);
213  if (not roGed.is_null())
214  xBvRoGed_ = rcp_dynamic_cast<BVROGED>(ged, true);
215  else if (not beLoc.is_null())
216  {
217  if (useTimeDerivativeSolutionVector_)
218  x_ = rcp_dynamic_cast<ProductVectorBase<double>>(beLoc->get_dxdt());
219  else // if (not useTimeDerivativeSolutionVector_)
220  x_ = rcp_dynamic_cast<ProductVectorBase<double>>(beLoc->get_x());
221  TEUCHOS_TEST_FOR_EXCEPTION(x_.is_null(), logic_error, "Gather " \
222  "Residual: Can't find the x_ vector in GEDC \"" << globalDataKey_ <<
223  "\" (" << post << "). A cast failed for " << ged << ". Type is " <<
224  typeName(*ged));
225  } // end if we have a roGed or beLoc
226  } // end of the old path
227 } // end of preEvaluate() (Residual Specialization)
228 
230 //
231 // evaluateFields() (Residual Specialization)
232 //
234 template<typename TRAITS, typename LO, typename GO>
235 void
236 panzer::
239  typename TRAITS::EvalData workset)
240 {
241  using PHX::MDField;
242  using std::size_t;
243  using std::string;
244  using std::vector;
245  using Teuchos::ArrayRCP;
246  using Teuchos::ptrFromRef;
247  using Teuchos::RCP;
248  using Teuchos::rcp_dynamic_cast;
250  using Thyra::SpmdVectorBase;
251  using Thyra::VectorBase;
252 
253  // For convenience, pull out some objects from the workset.
254  string blockId(this->wda(workset).block_id);
255  const vector<size_t>& localCellIds = this->wda(workset).cell_local_ids;
256  int numFields(gatherFields_.size()), numCells(localCellIds.size());
257 
258  // Loop over the fields to be gathered.
259  for (int fieldInd(0); fieldInd < numFields; ++fieldInd)
260  {
261  MDField<ScalarT, Cell, NODE>& field = gatherFields_[fieldInd];
262  auto field_h = Kokkos::create_mirror_view(field.get_static_view());
263 
264  int indexerId(indexerIds_[fieldInd]),
265  subFieldNum(subFieldIds_[fieldInd]);
266 
267  // Grab the local data for inputing.
268  ArrayRCP<const double> x;
269  Teuchos::RCP<const ReadOnlyVector_GlobalEvaluationData> xEvRoGed;
270 
271  if(not x_.is_null()) {
272  rcp_dynamic_cast<SpmdVectorBase<double>>(x_->
273  getNonconstVectorBlock(indexerId))->getLocalData(ptrFromRef(x));
274  }
275  else {
276  xEvRoGed = xBvRoGed_->getGEDBlock(indexerId);
277  }
278 
279  auto subRowIndexer = indexers_[indexerId];
280  const vector<int>& elmtOffset =
281  subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
282  int numBases(elmtOffset.size());
283 
284  auto LIDs = subRowIndexer->getLIDs();
285  auto LIDs_h = Kokkos::create_mirror_view(LIDs);
286  Kokkos::deep_copy(LIDs_h, LIDs);
287 
288  // Gather operation for each cell in the workset.
289  for (int cell(0); cell < numCells; ++cell)
290  {
291  LO cellLocalId = localCellIds[cell];
292 
293  // Loop over the basis functions and fill the fields.
294  for (int basis(0); basis < numBases; ++basis)
295  {
296  int offset(elmtOffset[basis]), lid(LIDs_h(cellLocalId, offset));
297  if(x_.is_null())
298  field_h(cell, basis) = (*xEvRoGed)[lid];
299  else
300  field_h(cell, basis) = x[lid];
301  } // end loop over the basis functions
302  } // end loop over localCellIds
303  Kokkos::deep_copy(field.get_static_view(), field_h);
304  } // end loop over the fields to be gathered
305 } // end of evaluateFields() (Residual Specialization)
306 
308 //
309 // Initializing Constructor (Tangent Specialization)
310 //
312 template<typename TRAITS, typename LO, typename GO>
315  const std::vector<Teuchos::RCP<const GlobalIndexer>>&
316  indexers,
317  const Teuchos::ParameterList& p)
318  :
319  indexers_(indexers),
320  hasTangentFields_(false)
321 {
322  using panzer::PureBasis;
323  using PHX::MDField;
324  using PHX::print;
325  using std::size_t;
326  using std::string;
327  using std::vector;
328  using Teuchos::RCP;
329  using vvstring = std::vector<std::vector<std::string>>;
330  GatherSolution_Input input;
331  input.setParameterList(p);
332  const vector<string>& names = input.getDofNames();
333  RCP<const PureBasis> basis = input.getBasis();
334  const vvstring& tangentFieldNames = input.getTangentNames();
335  indexerNames_ = input.getIndexerNames();
336  useTimeDerivativeSolutionVector_ = input.useTimeDerivativeSolutionVector();
337  globalDataKey_ = input.getGlobalDataKey();
338 
339  // Allocate the fields.
340  int numFields(names.size());
341  gatherFields_.resize(numFields);
342  for (int fd(0); fd < numFields; ++fd)
343  {
344  gatherFields_[fd] =
345  MDField<ScalarT, Cell, NODE>(names[fd], basis->functional);
346  this->addEvaluatedField(gatherFields_[fd]);
347  } // end loop over names
348 
349  // Set up the dependent tangent fields, if requested.
350  if (tangentFieldNames.size() > 0)
351  {
352  TEUCHOS_ASSERT(gatherFields_.size() == tangentFieldNames.size());
353  hasTangentFields_ = true;
354  tangentFields_.resize(numFields);
355  for (int fd(0); fd < numFields; ++fd)
356  {
357  int numTangentFields(tangentFieldNames[fd].size());
358  tangentFields_[fd].resize(numTangentFields);
359  for (int i(0); i < numTangentFields; ++i)
360  {
361  tangentFields_[fd][i] =
362  MDField<const ScalarT, Cell, NODE>(tangentFieldNames[fd][i],
363  basis->functional);
364  this->addDependentField(tangentFields_[fd][i]);
365  } // end loop over tangentFieldNames
366  } // end loop over gatherFields_
367  } // end if we have tangent fields
368 
369  // Figure out what the first active name is.
370  string firstName("<none>");
371  if (numFields > 0)
372  firstName = names[0];
373  string n("GatherSolution Tangent (BlockedEpetra): " + firstName + " (" +
374  print<EvalT>() + ")");
375  this->setName(n);
376 } // end of Initializing Constructor (Tangent Specialization)
377 
379 //
380 // postRegistrationSetup() (Tangent Specialization)
381 //
383 template<typename TRAITS, typename LO, typename GO>
384 void
387  typename TRAITS::SetupData /* d */,
388  PHX::FieldManager<TRAITS>& /* fm */)
389 {
390  using std::size_t;
391  using std::string;
392  TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size());
393  int numFields(gatherFields_.size());
394  indexerIds_.resize(numFields);
395  subFieldIds_.resize(numFields);
396  for (int fd(0); fd < numFields; ++fd)
397  {
398  // Get the field ID from the DOF manager.
399  const string& fieldName(indexerNames_[fd]);
400  indexerIds_[fd] = getFieldBlock(fieldName, indexers_);
401  subFieldIds_[fd] = indexers_[indexerIds_[fd]]->getFieldNum(fieldName);
402  TEUCHOS_ASSERT(indexerIds_[fd] >= 0);
403  } // end loop over gatherFields_
404  indexerNames_.clear();
405 } // end of postRegistrationSetup() (Tangent Specialization)
406 
408 //
409 // preEvaluate() (Tangent Specialization)
410 //
412 template<typename TRAITS, typename LO, typename GO>
413 void
416  typename TRAITS::PreEvalData d)
417 {
418  using std::logic_error;
419  using std::string;
420  using Teuchos::RCP;
421  using Teuchos::rcp_dynamic_cast;
422  using Teuchos::typeName;
426  using GED = panzer::GlobalEvaluationData;
427 
428  // First try the refactored ReadOnly container.
429  RCP<GED> ged;
430  string post(useTimeDerivativeSolutionVector_ ? " - Xdot" : " - X");
431  if (d.gedc->containsDataObject(globalDataKey_ + post))
432  {
433  ged = d.gedc->getDataObject(globalDataKey_ + post);
434  xBvRoGed_ = rcp_dynamic_cast<BVROGED>(ged, true);
435  return;
436  } // end of the refactored ReadOnly way
437 
438  // Now try the old path.
439  {
440  ged = d.gedc->getDataObject(globalDataKey_);
441 
442  // Extract the linear object container.
443  auto roGed = rcp_dynamic_cast<const BVROGED>(ged);
444  auto beLoc = rcp_dynamic_cast<const BELOC>(ged);
445  if (not roGed.is_null())
446  xBvRoGed_ = rcp_dynamic_cast<BVROGED>(ged, true);
447  else if (not beLoc.is_null())
448  {
449  if (useTimeDerivativeSolutionVector_)
450  x_ = rcp_dynamic_cast<ProductVectorBase<double>>(beLoc->get_dxdt());
451  else // if (not useTimeDerivativeSolutionVector_)
452  x_ = rcp_dynamic_cast<ProductVectorBase<double>>(beLoc->get_x());
453  TEUCHOS_TEST_FOR_EXCEPTION(x_.is_null(), logic_error, "Gather " \
454  "Tangent: Can't find the x_ vector in GEDC \"" << globalDataKey_ <<
455  "\" (" << post << "). A cast failed for " << ged << ". Type is " <<
456  typeName(*ged));
457  } // end if we have a roGed or beLoc
458  } // end of the old path
459 } // end of preEvaluate() (Tangent Specialization)
460 
462 //
463 // evaluateFields() (Tangent Specialization)
464 //
466 template<typename TRAITS, typename LO, typename GO>
467 void
470  typename TRAITS::EvalData workset)
471 {
472  using PHX::MDField;
473  using std::pair;
474  using std::size_t;
475  using std::string;
476  using std::vector;
477  using Teuchos::ArrayRCP;
478  using Teuchos::ptrFromRef;
479  using Teuchos::RCP;
480  using Teuchos::rcp_dynamic_cast;
482  using Thyra::SpmdVectorBase;
483  using Thyra::VectorBase;
484 
485  // For convenience, pull out some objects from the workset.
486  vector<pair<int, GO>> GIDs;
487  string blockId(this->wda(workset).block_id);
488  const vector<size_t>& localCellIds = this->wda(workset).cell_local_ids;
489  int numFields(gatherFields_.size()), numCells(localCellIds.size());
490 
491  if (x_.is_null())
492  {
493  // Loop over the fields to be gathered.
494  for (int fieldInd(0); fieldInd < numFields; ++fieldInd)
495  {
496  MDField<ScalarT, Cell, NODE>& field = gatherFields_[fieldInd];
497  int indexerId(indexerIds_[fieldInd]),
498  subFieldNum(subFieldIds_[fieldInd]);
499 
500  // Grab the local data for inputing.
501  auto xEvRoGed = xBvRoGed_->getGEDBlock(indexerId);
502  auto subRowIndexer = indexers_[indexerId];
503  const vector<int>& elmtOffset =
504  subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
505  int numBases(elmtOffset.size());
506 
507  // Gather operation for each cell in the workset.
508  for (int cell(0); cell < numCells; ++cell)
509  {
510  LO cellLocalId = localCellIds[cell];
511  auto LIDs = subRowIndexer->getElementLIDs(cellLocalId);
512 
513  // Loop over the basis functions and fill the fields.
514  for (int basis(0); basis < numBases; ++basis)
515  {
516  int offset(elmtOffset[basis]), lid(LIDs[offset]);
517  field(cell, basis) = (*xEvRoGed)[lid];
518  } // end loop over the basis functions
519  } // end loop over localCellIds
520  } // end loop over the fields to be gathered
521  }
522  else // if (not x_.is_null())
523  {
524  // Loop over the fields to be gathered.
525  for (int fieldInd(0); fieldInd < numFields; ++fieldInd)
526  {
527  MDField<ScalarT, Cell, NODE>& field = gatherFields_[fieldInd];
528  int indexerId(indexerIds_[fieldInd]),
529  subFieldNum(subFieldIds_[fieldInd]);
530 
531  // Grab the local data for inputing.
532  ArrayRCP<const double> x;
533  rcp_dynamic_cast<SpmdVectorBase<double>>(x_->
534  getNonconstVectorBlock(indexerId))->getLocalData(ptrFromRef(x));
535  auto subRowIndexer = indexers_[indexerId];
536  const vector<int>& elmtOffset =
537  subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
538  int numBases(elmtOffset.size());
539 
540  // Gather operation for each cell in the workset.
541  for (int cell(0); cell < numCells; ++cell)
542  {
543  LO cellLocalId = localCellIds[cell];
544  auto LIDs = subRowIndexer->getElementLIDs(cellLocalId);
545 
546  // Loop over the basis functions and fill the fields.
547  for (int basis(0); basis < numBases; ++basis)
548  {
549  int offset(elmtOffset[basis]), lid(LIDs[offset]);
550  field(cell, basis) = x[lid];
551  } // end loop over the basis functions
552  } // end loop over localCellIds
553  } // end loop over the fields to be gathered
554  } // end if (x_.is_null()) or not
555 
556  // Deal with the tangent fields.
557  if (hasTangentFields_)
558  {
559  // Loop over the fields to be gathered.
560  for (int fieldInd(0); fieldInd < numFields; ++fieldInd)
561  {
562  MDField<ScalarT, Cell, NODE>& field = gatherFields_[fieldInd];
563  int indexerId(indexerIds_[fieldInd]),
564  subFieldNum(subFieldIds_[fieldInd]);
565  auto subRowIndexer = indexers_[indexerId];
566  const vector<int>& elmtOffset =
567  subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
568  int numBases(elmtOffset.size());
569 
570  // Gather operation for each cell in the workset.
571  for (int cell(0); cell < numCells; ++cell)
572  {
573  // Loop over the basis functions and fill the fields.
574  for (int basis(0); basis < numBases; ++basis)
575  {
576  int numTangentFields(tangentFields_[fieldInd].size());
577  for (int i(0); i < numTangentFields; ++i)
578  field(cell, basis).fastAccessDx(i) =
579  tangentFields_[fieldInd][i](cell, basis).val();
580  } // end loop over the basis functions
581  } // end loop over localCellIds
582  } // end loop over the fields to be gathered
583  } // end if (hasTangentFields_)
584 } // end of evaluateFields() (Tangent Specialization)
585 
587 //
588 // Initializing Constructor (Jacobian Specialization)
589 //
591 template<typename TRAITS, typename LO, typename GO>
592 panzer::
595  const std::vector<Teuchos::RCP<const GlobalIndexer>>&
596  indexers,
597  const Teuchos::ParameterList& p)
598  :
599  indexers_(indexers)
600 {
601  using panzer::PureBasis;
602  using PHX::MDField;
603  using PHX::print;
604  using std::size_t;
605  using std::string;
606  using std::vector;
607  using Teuchos::RCP;
608  GatherSolution_Input input;
609  input.setParameterList(p);
610  const vector<string>& names = input.getDofNames();
611  RCP<const PureBasis> basis = input.getBasis();
612  indexerNames_ = input.getIndexerNames();
613  useTimeDerivativeSolutionVector_ = input.useTimeDerivativeSolutionVector();
614  globalDataKey_ = input.getGlobalDataKey();
615  gatherSeedIndex_ = input.getGatherSeedIndex();
616  sensitivitiesName_ = input.getSensitivitiesName();
617  disableSensitivities_ = not input.firstSensitivitiesAvailable();
618 
619  // Allocate the fields.
620  int numFields(names.size());
621  gatherFields_.resize(numFields);
622  for (int fd(0); fd < numFields; ++fd)
623  {
624  MDField<ScalarT, Cell, NODE> f(names[fd], basis->functional);
625  gatherFields_[fd] = f;
626  this->addEvaluatedField(gatherFields_[fd]);
627  } // end loop over names
628 
629  // Figure out what the first active name is.
630  string firstName("<none>"), n("GatherSolution (BlockedEpetra");
631  if (numFields > 0)
632  firstName = names[0];
633  if (disableSensitivities_)
634  n += ", No Sensitivities";
635  n += "): " + firstName + " (" + print<EvalT>() + ")";
636  this->setName(n);
637 } // end of Initializing Constructor (Jacobian Specialization)
638 
640 //
641 // postRegistrationSetup() (Jacobian Specialization)
642 //
644 template<typename TRAITS, typename LO, typename GO>
645 void
646 panzer::
649  typename TRAITS::SetupData /* d */,
650  PHX::FieldManager<TRAITS>& /* fm */)
651 {
652  using std::size_t;
653  using std::string;
654  TEUCHOS_ASSERT(gatherFields_.size() == indexerNames_.size());
655  int numFields(gatherFields_.size());
656  indexerIds_.resize(numFields);
657  subFieldIds_.resize(numFields);
658  for (int fd(0); fd < numFields; ++fd)
659  {
660  // Get the field ID from the DOF manager.
661  const string& fieldName(indexerNames_[fd]);
662  indexerIds_[fd] = getFieldBlock(fieldName, indexers_);
663  subFieldIds_[fd] = indexers_[indexerIds_[fd]]->getFieldNum(fieldName);
664  TEUCHOS_ASSERT(indexerIds_[fd] >= 0);
665  } // end loop over gatherFields_
666  indexerNames_.clear();
667 } // end of postRegistrationSetup() (Jacobian Specialization)
668 
670 //
671 // preEvaluate() (Jacobian Specialization)
672 //
674 template<typename TRAITS, typename LO, typename GO>
675 void
676 panzer::
679  typename TRAITS::PreEvalData d)
680 {
681  using std::endl;
682  using std::logic_error;
683  using std::string;
684  using Teuchos::RCP;
685  using Teuchos::rcp_dynamic_cast;
686  using Teuchos::typeName;
690  using GED = panzer::GlobalEvaluationData;
691 
692  // Manage sensitivities.
693  applySensitivities_ = false;
694  if ((not disableSensitivities_ ) and
695  (d.first_sensitivities_name == sensitivitiesName_))
696  applySensitivities_ = true;
697 
698  // First try the refactored ReadOnly container.
699  RCP<GED> ged;
700  string post(useTimeDerivativeSolutionVector_ ? " - Xdot" : " - X");
701  if (d.gedc->containsDataObject(globalDataKey_ + post))
702  {
703  ged = d.gedc->getDataObject(globalDataKey_ + post);
704  xBvRoGed_ = rcp_dynamic_cast<BVROGED>(ged, true);
705  return;
706  } // end of the refactored ReadOnly way
707 
708  // Now try the old path.
709  {
710  ged = d.gedc->getDataObject(globalDataKey_);
711 
712  // Extract the linear object container.
713  auto roGed = rcp_dynamic_cast<const BVROGED>(ged);
714  auto beLoc = rcp_dynamic_cast<const BELOC>(ged);
715  if (not roGed.is_null())
716  xBvRoGed_ = rcp_dynamic_cast<BVROGED>(ged, true);
717  else if (not beLoc.is_null())
718  {
719  if (useTimeDerivativeSolutionVector_)
720  x_ = rcp_dynamic_cast<ProductVectorBase<double>>(beLoc->get_dxdt());
721  else // if (not useTimeDerivativeSolutionVector_)
722  x_ = rcp_dynamic_cast<ProductVectorBase<double>>(beLoc->get_x());
723  TEUCHOS_TEST_FOR_EXCEPTION(x_.is_null(), logic_error, "Gather " \
724  "Jacobian: Can't find x vector in GEDC \"" << globalDataKey_ <<
725  "\" (" << post << "). A cast failed for " << ged << ". Type is " <<
726  typeName(*ged));
727  } // end if we have a roGed or beLoc
728  } // end of the old path
729 } // end of preEvaluate() (Jacobian Specialization)
730 
732 //
733 // evaluateFields() (Jacobian Specialization)
734 //
736 template<typename TRAITS, typename LO, typename GO>
737 void
738 panzer::
741  typename TRAITS::EvalData workset)
742 {
743  using PHX::MDField;
744  using std::size_t;
745  using std::string;
746  using std::vector;
747  using Teuchos::ArrayRCP;
748  using Teuchos::ptrFromRef;
749  using Teuchos::RCP;
750  using Teuchos::rcp_dynamic_cast;
752  using Thyra::SpmdVectorBase;
753  using Thyra::VectorBase;
754 
755  // For convenience, pull out some objects from the workset.
756  string blockId(this->wda(workset).block_id);
757  const vector<size_t>& localCellIds = this->wda(workset).cell_local_ids;
758  int numFields(gatherFields_.size()), numCells(localCellIds.size());
759  double seedValue(0);
760  if (applySensitivities_)
761  {
762  if ((useTimeDerivativeSolutionVector_) and (gatherSeedIndex_ < 0))
763  seedValue = workset.alpha;
764  else if (gatherSeedIndex_ < 0)
765  seedValue = workset.beta;
766  else if (not useTimeDerivativeSolutionVector_)
767  seedValue = workset.gather_seeds[gatherSeedIndex_];
768  else
769  TEUCHOS_ASSERT(false);
770  } // end if (applySensitivities_)
771 
772  // Turn off sensitivies: This may be faster if we don't expand the term, but
773  // I suspect not, because anywhere it is used the full complement of
774  // sensitivies will be needed anyway.
775  if (not applySensitivities_)
776  seedValue = 0.0;
777 
778  vector<int> blockOffsets;
779  computeBlockOffsets(blockId, indexers_, blockOffsets);
780  int numDerivs(blockOffsets[blockOffsets.size() - 1]);
781 
782  // NOTE: A reordering of these loops will likely improve performance. The
783  // "getGIDFieldOffsets may be expensive. However the "getElementGIDs"
784  // can be cheaper. However the lookup for LIDs may be more expensive!
785 
786  // Loop over the fields to be gathered.
787  for (int fieldInd(0); fieldInd < numFields; ++fieldInd)
788  {
789  MDField<ScalarT, Cell, NODE>& field = gatherFields_[fieldInd];
790  auto field_h = Kokkos::create_mirror_view(field.get_view());
791 
792  int indexerId(indexerIds_[fieldInd]), subFieldNum(subFieldIds_[fieldInd]);
793 
794  // Grab the local data for inputing.
795  ArrayRCP<const double> x;
796  Teuchos::RCP<const ReadOnlyVector_GlobalEvaluationData> xEvRoGed;
797  if(not x_.is_null()) {
798  rcp_dynamic_cast<SpmdVectorBase<double>>(x_->getNonconstVectorBlock(indexerId))->getLocalData(ptrFromRef(x));
799  }else {
800  xEvRoGed = xBvRoGed_->getGEDBlock(indexerId);
801  }
802 
803  auto subRowIndexer = indexers_[indexerId];
804  const vector<int>& elmtOffset =
805  subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
806  int startBlkOffset(blockOffsets[indexerId]), numBases(elmtOffset.size());
807 
808  auto LIDs = subRowIndexer->getLIDs();
809  auto LIDs_h = Kokkos::create_mirror_view(LIDs);
810  Kokkos::deep_copy(LIDs_h, LIDs);
811 
812  // Gather operation for each cell in the workset.
813  for (int cell(0); cell < numCells; ++cell)
814  {
815  LO cellLocalId = localCellIds[cell];
816 
817  // Loop over the basis functions and fill the fields.
818  for (int basis(0); basis < numBases; ++basis)
819  {
820  // Set the value and seed the FAD object.
821  int offset(elmtOffset[basis]), lid(LIDs_h(cellLocalId, offset));
822  if(x_.is_null())
823  field_h(cell, basis) = ScalarT(numDerivs, (*xEvRoGed)[lid]);
824  else
825  field_h(cell, basis) = ScalarT(numDerivs, x[lid]);
826 
827  field_h(cell, basis).fastAccessDx(startBlkOffset + offset) = seedValue;
828  } // end loop over the basis functions
829  } // end loop over localCellIds
830  Kokkos::deep_copy(field.get_static_view(), field_h);
831  } // end loop over the fields to be gathered
832 } // end of evaluateFields() (Jacobian Specialization)
833 
834 #endif // __Panzer_GatherSolution_BlockedEpetra_impl_hpp__
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)
void setParameterList(const Teuchos::ParameterList &pl)
const std::vector< std::string > & getDofNames() const
The names of the DOFs to be gathered (all types)
void computeBlockOffsets(const std::string &blockId, const std::vector< Teuchos::RCP< GlobalIndexer >> &ugis, std::vector< int > &blockOffsets)
int numFields
This class encapsulates the needs of a gather operation to do a halo exchange for blocked vectors...
int getGatherSeedIndex() const
What index to use for initializing the seed (Jacobian and Hessian)
bool useTimeDerivativeSolutionVector() const
Gather a time derivative vector? (all types)
int getFieldBlock(const std::string &fieldName, const std::vector< Teuchos::RCP< const GlobalIndexer >> &ugis)
bool firstSensitivitiesAvailable()
Are first derivative sensitivities enabled or disabled? (Jacobian and Hessian)
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &fm)
Post-Registration Setup.
GatherSolution_BlockedEpetra(const Teuchos::ParameterList &p)
Constructor.
std::string getSensitivitiesName() const
The name of the sensitivities. Enables sensitivities at "preEvaluate" time (Jacobian and Hessian) ...
Gathers solution values from the Newton solution vector into the nodal fields of the field manager...
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
void evaluateFields(typename TRAITS::EvalData d)
Evaluate Fields.
Description and data layouts associated with a particular basis.
const std::vector< std::vector< std::string > > & getTangentNames() const
Get the name of the tangent fields (tangent only)