Panzer  Version of the Day
Panzer_ScatterResidual_Tpetra_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_SCATTER_RESIDUAL_TPETRA_IMPL_HPP
44 #define PANZER_SCATTER_RESIDUAL_TPETRA_IMPL_HPP
45 
46 #include "Teuchos_RCP.hpp"
47 #include "Teuchos_Assert.hpp"
48 
49 #include "Phalanx_DataLayout.hpp"
50 
51 #include "Panzer_GlobalIndexer.hpp"
52 #include "Panzer_PureBasis.hpp"
57 
58 #include "Phalanx_DataLayout_MDALayout.hpp"
59 
60 #include "Teuchos_FancyOStream.hpp"
61 
62 #include "Tpetra_Vector.hpp"
63 #include "Tpetra_Map.hpp"
64 #include "Tpetra_CrsMatrix.hpp"
65 
66 // **********************************************************************
67 // Specialization: Residual
68 // **********************************************************************
69 
70 template<typename TRAITS,typename LO,typename GO,typename NodeT>
72 ScatterResidual_Tpetra(const Teuchos::RCP<const panzer::GlobalIndexer> & indexer,
73  const Teuchos::ParameterList& p)
74  : globalIndexer_(indexer)
75  , globalDataKey_("Residual Scatter Container")
76 {
77  std::string scatterName = p.get<std::string>("Scatter Name");
78  scatterHolder_ =
79  Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
80 
81  // get names to be evaluated
82  const std::vector<std::string>& names =
83  *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
84 
85  // grab map from evaluated names to field names
86  fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
87 
88  Teuchos::RCP<PHX::DataLayout> dl =
89  p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional;
90 
91  // build the vector of fields that this is dependent on
92  scatterFields_.resize(names.size());
93  scratch_offsets_.resize(names.size());
94  for (std::size_t eq = 0; eq < names.size(); ++eq) {
95  scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
96 
97  // tell the field manager that we depend on this field
98  this->addDependentField(scatterFields_[eq]);
99  }
100 
101  // this is what this evaluator provides
102  this->addEvaluatedField(*scatterHolder_);
103 
104  if (p.isType<std::string>("Global Data Key"))
105  globalDataKey_ = p.get<std::string>("Global Data Key");
106 
107  this->setName(scatterName+" Scatter Residual");
108 }
109 
110 // **********************************************************************
111 template<typename TRAITS,typename LO,typename GO,typename NodeT>
113 postRegistrationSetup(typename TRAITS::SetupData d,
114  PHX::FieldManager<TRAITS>& /* fm */)
115 {
116  fieldIds_.resize(scatterFields_.size());
117  const Workset & workset_0 = (*d.worksets_)[0];
118  std::string blockId = this->wda(workset_0).block_id;
119 
120 
121  // load required field numbers for fast use
122  for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
123  // get field ID from DOF manager
124  std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
125  fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
126 
127  const std::vector<int> & offsets = globalIndexer_->getGIDFieldOffsets(blockId,fieldIds_[fd]);
128  scratch_offsets_[fd] = PHX::View<int*>("offsets",offsets.size());
129  Kokkos::deep_copy(scratch_offsets_[fd], Kokkos::View<const int*, Kokkos::HostSpace, Kokkos::MemoryUnmanaged>(offsets.data(), offsets.size()));
130  }
131  scratch_lids_ = PHX::View<LO**>("lids",scatterFields_[0].extent(0),
132  globalIndexer_->getElementBlockGIDCount(blockId));
133 
134 }
135 
136 // **********************************************************************
137 template<typename TRAITS,typename LO,typename GO,typename NodeT>
139 preEvaluate(typename TRAITS::PreEvalData d)
140 {
142 
143  // extract linear object container
144  tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(d.gedc->getDataObject(globalDataKey_));
145 
146  if(tpetraContainer_==Teuchos::null) {
147  // extract linear object container
148  Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),true)->getGhostedLOC();
149  tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(loc);
150  }
151 }
152 
153 
154 // **********************************************************************
155 // Specialization: Tangent
156 // **********************************************************************
157 
158 template<typename TRAITS,typename LO,typename GO,typename NodeT>
160 ScatterResidual_Tpetra(const Teuchos::RCP<const panzer::GlobalIndexer> & indexer,
161  const Teuchos::ParameterList& p)
162  : globalIndexer_(indexer)
163  , globalDataKey_("Residual Scatter Container")
164 {
165  std::string scatterName = p.get<std::string>("Scatter Name");
166  scatterHolder_ =
167  Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
168 
169  // get names to be evaluated
170  const std::vector<std::string>& names =
171  *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
172 
173  // grab map from evaluated names to field names
174  fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
175 
176  Teuchos::RCP<PHX::DataLayout> dl =
177  p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional;
178 
179  // build the vector of fields that this is dependent on
180  scatterFields_.resize(names.size());
181  for (std::size_t eq = 0; eq < names.size(); ++eq) {
182  scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
183 
184  // tell the field manager that we depend on this field
185  this->addDependentField(scatterFields_[eq]);
186  }
187 
188  // this is what this evaluator provides
189  this->addEvaluatedField(*scatterHolder_);
190 
191  if (p.isType<std::string>("Global Data Key"))
192  globalDataKey_ = p.get<std::string>("Global Data Key");
193 
194  this->setName(scatterName+" Scatter Tangent");
195 }
196 
197 // **********************************************************************
198 template<typename TRAITS,typename LO,typename GO,typename NodeT>
200 postRegistrationSetup(typename TRAITS::SetupData /* d */,
201  PHX::FieldManager<TRAITS>& /* fm */)
202 {
203  fieldIds_.resize(scatterFields_.size());
204  // load required field numbers for fast use
205  for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
206  // get field ID from DOF manager
207  std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
208  fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
209  }
210 }
211 
212 // **********************************************************************
213 template<typename TRAITS,typename LO,typename GO,typename NodeT>
215 preEvaluate(typename TRAITS::PreEvalData d)
216 {
217  using Teuchos::RCP;
218  using Teuchos::rcp_dynamic_cast;
219 
221 
222  // this is the list of parameters and their names that this scatter has to account for
223  std::vector<std::string> activeParameters =
224  rcp_dynamic_cast<ParameterList_GlobalEvaluationData>(d.gedc->getDataObject("PARAMETER_NAMES"))->getActiveParameters();
225 
226  dfdp_vectors_.clear();
227  for(std::size_t i=0;i<activeParameters.size();i++) {
228  RCP<typename LOC::VectorType> vec =
229  rcp_dynamic_cast<LOC>(d.gedc->getDataObject(activeParameters[i]),true)->get_f();
230  Teuchos::ArrayRCP<double> vec_array = vec->get1dViewNonConst();
231  dfdp_vectors_.push_back(vec_array);
232  }
233 }
234 
235 // **********************************************************************
236 template<typename TRAITS,typename LO,typename GO,typename NodeT>
238 evaluateFields(typename TRAITS::EvalData workset)
239 {
240  // for convenience pull out some objects from workset
241  std::string blockId = this->wda(workset).block_id;
242  const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
243 
244  // NOTE: A reordering of these loops will likely improve performance
245  // The "getGIDFieldOffsets may be expensive. However the
246  // "getElementGIDs" can be cheaper. However the lookup for LIDs
247  // may be more expensive!
248 
249  // scatter operation for each cell in workset
250  for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
251  std::size_t cellLocalId = localCellIds[worksetCellIndex];
252 
253  auto LIDs = globalIndexer_->getElementLIDs(cellLocalId);
254 
255  // loop over each field to be scattered
256  for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
257  int fieldNum = fieldIds_[fieldIndex];
258  const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
259 
260  // loop over basis functions
261  for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
262  int offset = elmtOffset[basis];
263  LO lid = LIDs[offset];
264  ScalarT value = (scatterFields_[fieldIndex])(worksetCellIndex,basis);
265  for(int d=0;d<value.size();d++)
266  dfdp_vectors_[d][lid] += value.fastAccessDx(d);
267  }
268  }
269  }
270 }
271 
272 // **********************************************************************
273 // Specialization: Jacobian
274 // **********************************************************************
275 
276 template<typename TRAITS,typename LO,typename GO,typename NodeT>
278 ScatterResidual_Tpetra(const Teuchos::RCP<const GlobalIndexer> & indexer,
279  const Teuchos::ParameterList& p)
280  : globalIndexer_(indexer)
281  , globalDataKey_("Residual Scatter Container")
282  , my_derivative_size_(0)
283  , other_derivative_size_(0)
284 {
285  std::string scatterName = p.get<std::string>("Scatter Name");
286  scatterHolder_ =
287  Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
288 
289  // get names to be evaluated
290  const std::vector<std::string>& names =
291  *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
292 
293  // grab map from evaluated names to field names
294  fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
295 
296  Teuchos::RCP<PHX::DataLayout> dl =
297  p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional;
298 
299  // build the vector of fields that this is dependent on
300  scatterFields_.resize(names.size());
301  scratch_offsets_.resize(names.size());
302  for (std::size_t eq = 0; eq < names.size(); ++eq) {
303  scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
304 
305  // tell the field manager that we depend on this field
306  this->addDependentField(scatterFields_[eq]);
307  }
308 
309  // this is what this evaluator provides
310  this->addEvaluatedField(*scatterHolder_);
311 
312  if (p.isType<std::string>("Global Data Key"))
313  globalDataKey_ = p.get<std::string>("Global Data Key");
314 
315  this->setName(scatterName+" Scatter Residual (Jacobian)");
316 }
317 
318 // **********************************************************************
319 template<typename TRAITS,typename LO,typename GO,typename NodeT>
321 postRegistrationSetup(typename TRAITS::SetupData d,
322  PHX::FieldManager<TRAITS>& /* fm */)
323 {
324  fieldIds_.resize(scatterFields_.size());
325 
326  const Workset & workset_0 = (*d.worksets_)[0];
327  std::string blockId = this->wda(workset_0).block_id;
328 
329  // load required field numbers for fast use
330  for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
331  // get field ID from DOF manager
332  std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
333  fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
334 
335  int fieldNum = fieldIds_[fd];
336  const std::vector<int> & offsets = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
337  scratch_offsets_[fd] = PHX::View<int*>("offsets",offsets.size());
338  Kokkos::deep_copy(scratch_offsets_[fd], Kokkos::View<const int*, Kokkos::HostSpace, Kokkos::MemoryUnmanaged>(offsets.data(), offsets.size()));
339  }
340 
341  my_derivative_size_ = globalIndexer_->getElementBlockGIDCount(blockId);
342  if (Teuchos::nonnull(workset_0.other)) {
343  auto otherBlockId = workset_0.other->block_id;
344  other_derivative_size_ = globalIndexer_->getElementBlockGIDCount(otherBlockId);
345  }
346  scratch_lids_ = PHX::View<LO**>("lids",scatterFields_[0].extent(0),
347  my_derivative_size_ + other_derivative_size_);
348 }
349 
350 // **********************************************************************
351 template<typename TRAITS,typename LO,typename GO,typename NodeT>
353 preEvaluate(typename TRAITS::PreEvalData d)
354 {
356 
357  // extract linear object container
358  tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(d.gedc->getDataObject(globalDataKey_));
359 
360  if(tpetraContainer_==Teuchos::null) {
361  // extract linear object container
362  Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),true)->getGhostedLOC();
363  tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(loc);
364  }
365 }
366 
367 
368 // **********************************************************************
369 namespace panzer {
370 namespace {
371 
372 template <typename ScalarT,typename LO,typename GO,typename NodeT,typename LocalMatrixT>
373 class ScatterResidual_Jacobian_Functor {
374 public:
375  typedef typename PHX::Device execution_space;
376  typedef PHX::MDField<const ScalarT,Cell,NODE> FieldType;
377 
379  Kokkos::View<double**, Kokkos::LayoutLeft,PHX::Device> r_data;
380  LocalMatrixT jac; // Kokkos jacobian type
381 
382  PHX::View<const LO**> lids; // local indices for unknowns
383  PHX::View<const int*> offsets; // how to get a particular field
385 
386 
387  KOKKOS_INLINE_FUNCTION
388  void operator()(const unsigned int cell) const
389  {
390  LO cLIDs[256];
391  typename Sacado::ScalarType<ScalarT>::type vals[256];
392  int numIds = lids.extent(1);
393 
394  for(int i=0;i<numIds;i++)
395  cLIDs[i] = lids(cell,i);
396 
397  // loop over the basis functions (currently they are nodes)
398  for(std::size_t basis=0; basis < offsets.extent(0); basis++) {
399  typename FieldType::array_type::reference_type scatterField = field(cell,basis);
400  int offset = offsets(basis);
401  LO lid = lids(cell,offset);
402 
403  // Sum residual
404  if(fillResidual)
405  Kokkos::atomic_add(&r_data(lid,0), scatterField.val());
406 
407  // loop over the sensitivity indices: all DOFs on a cell
408  for(int sensIndex=0;sensIndex<numIds;++sensIndex)
409  vals[sensIndex] = scatterField.fastAccessDx(sensIndex);
410 
411  // Sum Jacobian
412  jac.sumIntoValues(lid, cLIDs,numIds, vals, true, true);
413  } // end basis
414  }
415 };
416 
417 template <typename ScalarT,typename LO,typename GO,typename NodeT>
418 class ScatterResidual_Residual_Functor {
419 public:
420  typedef typename PHX::Device execution_space;
421  typedef PHX::MDField<const ScalarT,Cell,NODE> FieldType;
422 
423  Kokkos::View<double**, Kokkos::LayoutLeft,PHX::Device> r_data;
424 
425  PHX::View<const LO**> lids; // local indices for unknowns
426  PHX::View<const int*> offsets; // how to get a particular field
428 
429 
430  KOKKOS_INLINE_FUNCTION
431  void operator()(const unsigned int cell) const
432  {
433 
434  // loop over the basis functions (currently they are nodes)
435  for(std::size_t basis=0; basis < offsets.extent(0); basis++) {
436  int offset = offsets(basis);
437  LO lid = lids(cell,offset);
438  Kokkos::atomic_add(&r_data(lid,0), field(cell,basis));
439 
440  } // end basis
441  }
442 };
443 
444 }
445 }
446 
447 // **********************************************************************
448 template<typename TRAITS,typename LO,typename GO,typename NodeT>
450 evaluateFields(typename TRAITS::EvalData workset)
451 {
453 
454  // for convenience pull out some objects from workset
455  std::string blockId = this->wda(workset).block_id;
456 
457  Teuchos::RCP<typename LOC::VectorType> r = tpetraContainer_->get_f();
458 
459  globalIndexer_->getElementLIDs(this->wda(workset).cell_local_ids_k,scratch_lids_);
460 
461  ScatterResidual_Residual_Functor<ScalarT,LO,GO,NodeT> functor;
462  functor.r_data = r->getLocalViewDevice(Tpetra::Access::ReadWrite);
463  functor.lids = scratch_lids_;
464 
465  // for each field, do a parallel for loop
466  for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
467  functor.offsets = scratch_offsets_[fieldIndex];
468  functor.field = scatterFields_[fieldIndex];
469 
470  Kokkos::parallel_for(workset.num_cells,functor);
471  }
472 }
473 
474 // **********************************************************************
475 template<typename TRAITS,typename LO,typename GO,typename NodeT>
477 evaluateFields(typename TRAITS::EvalData workset)
478 {
480 
481  typedef typename LOC::CrsMatrixType::local_matrix_device_type LocalMatrixT;
482 
483  // for convenience pull out some objects from workset
484  std::string blockId = this->wda(workset).block_id;
485 
486  Teuchos::RCP<typename LOC::VectorType> r = tpetraContainer_->get_f();
487  Teuchos::RCP<typename LOC::CrsMatrixType> Jac = tpetraContainer_->get_A();
488 
489  // Cache scratch lids. For interface bc problems the derivative
490  // dimension extent spans two cells. Use subviews to get the self
491  // lids and the other lids.
492  if (Teuchos::nonnull(workset.other)) {
493  auto my_scratch_lids = Kokkos::subview(scratch_lids_,Kokkos::ALL,std::make_pair(0,my_derivative_size_));
494  globalIndexer_->getElementLIDs(this->wda(workset).cell_local_ids_k,my_scratch_lids);
495  auto other_scratch_lids = Kokkos::subview(scratch_lids_,Kokkos::ALL,std::make_pair(my_derivative_size_,my_derivative_size_ + other_derivative_size_));
496  globalIndexer_->getElementLIDs(workset.other->cell_local_ids_k,other_scratch_lids);
497  }
498  else {
499  globalIndexer_->getElementLIDs(this->wda(workset).cell_local_ids_k,scratch_lids_);
500  }
501 
502  ScatterResidual_Jacobian_Functor<ScalarT,LO,GO,NodeT,LocalMatrixT> functor;
503  functor.fillResidual = (r!=Teuchos::null);
504  if(functor.fillResidual)
505  functor.r_data = r->getLocalViewDevice(Tpetra::Access::ReadWrite);
506  functor.jac = Jac->getLocalMatrixDevice();
507  functor.lids = scratch_lids_;
508 
509  // for each field, do a parallel for loop
510  for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
511  functor.offsets = scratch_offsets_[fieldIndex];
512  functor.field = scatterFields_[fieldIndex];
513 
514  Kokkos::parallel_for(workset.num_cells,functor);
515  }
516 
517 }
518 
519 // **********************************************************************
520 
521 #endif
FieldType
The type of discretization to use for a field pattern.
Teuchos::RCP< WorksetDetails > other
Kokkos::View< double **, Kokkos::LayoutLeft, PHX::Device > r_data
std::string block_id
DEPRECATED - use: getElementBlock()
Pushes residual values into the residual vector for a Newton-based solve.
PHX::View< const int * > offsets
const Teuchos::RCP< VectorType > get_f() const
PHX::View< const LO ** > lids