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