Panzer  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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"
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>
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 
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] = Kokkos::View<int*,PHX::Device>("offsets",offsets.size());
129  for(std::size_t i=0;i<offsets.size();i++)
130  scratch_offsets_[fd](i) = offsets[i];
131  }
132  scratch_lids_ = Kokkos::View<LO**,PHX::Device>("lids",scatterFields_[0].extent(0),
133  globalIndexer_->getElementBlockGIDCount(blockId));
134 
135 }
136 
137 // **********************************************************************
138 template<typename TRAITS,typename LO,typename GO,typename NodeT>
140 preEvaluate(typename TRAITS::PreEvalData d)
141 {
143 
144  // extract linear object container
145  tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(d.gedc->getDataObject(globalDataKey_));
146 
147  if(tpetraContainer_==Teuchos::null) {
148  // extract linear object container
149  Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),true)->getGhostedLOC();
150  tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(loc);
151  }
152 }
153 
154 
155 // **********************************************************************
156 // Specialization: Tangent
157 // **********************************************************************
158 
159 template<typename TRAITS,typename LO,typename GO,typename NodeT>
162  const Teuchos::ParameterList& p)
163  : globalIndexer_(indexer)
164  , globalDataKey_("Residual Scatter Container")
165 {
166  std::string scatterName = p.get<std::string>("Scatter Name");
167  scatterHolder_ =
168  Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
169 
170  // get names to be evaluated
171  const std::vector<std::string>& names =
172  *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
173 
174  // grab map from evaluated names to field names
175  fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
176 
178  p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional;
179 
180  // build the vector of fields that this is dependent on
181  scatterFields_.resize(names.size());
182  for (std::size_t eq = 0; eq < names.size(); ++eq) {
183  scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
184 
185  // tell the field manager that we depend on this field
186  this->addDependentField(scatterFields_[eq]);
187  }
188 
189  // this is what this evaluator provides
190  this->addEvaluatedField(*scatterHolder_);
191 
192  if (p.isType<std::string>("Global Data Key"))
193  globalDataKey_ = p.get<std::string>("Global Data Key");
194 
195  this->setName(scatterName+" Scatter Tangent");
196 }
197 
198 // **********************************************************************
199 template<typename TRAITS,typename LO,typename GO,typename NodeT>
201 postRegistrationSetup(typename TRAITS::SetupData /* d */,
202  PHX::FieldManager<TRAITS>& /* fm */)
203 {
204  fieldIds_.resize(scatterFields_.size());
205  // load required field numbers for fast use
206  for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
207  // get field ID from DOF manager
208  std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
209  fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
210  }
211 }
212 
213 // **********************************************************************
214 template<typename TRAITS,typename LO,typename GO,typename NodeT>
216 preEvaluate(typename TRAITS::PreEvalData d)
217 {
218  using Teuchos::RCP;
219  using Teuchos::rcp_dynamic_cast;
220 
222 
223  // this is the list of parameters and their names that this scatter has to account for
224  std::vector<std::string> activeParameters =
225  rcp_dynamic_cast<ParameterList_GlobalEvaluationData>(d.gedc->getDataObject("PARAMETER_NAMES"))->getActiveParameters();
226 
227  dfdp_vectors_.clear();
228  for(std::size_t i=0;i<activeParameters.size();i++) {
230  rcp_dynamic_cast<LOC>(d.gedc->getDataObject(activeParameters[i]),true)->get_f();
231  Teuchos::ArrayRCP<double> vec_array = vec->get1dViewNonConst();
232  dfdp_vectors_.push_back(vec_array);
233  }
234 }
235 
236 // **********************************************************************
237 template<typename TRAITS,typename LO,typename GO,typename NodeT>
239 evaluateFields(typename TRAITS::EvalData workset)
240 {
241  // for convenience pull out some objects from workset
242  std::string blockId = this->wda(workset).block_id;
243  const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
244 
245  // NOTE: A reordering of these loops will likely improve performance
246  // The "getGIDFieldOffsets may be expensive. However the
247  // "getElementGIDs" can be cheaper. However the lookup for LIDs
248  // may be more expensive!
249 
250  // scatter operation for each cell in workset
251  for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
252  std::size_t cellLocalId = localCellIds[worksetCellIndex];
253 
254  auto LIDs = globalIndexer_->getElementLIDs(cellLocalId);
255 
256  // loop over each field to be scattered
257  for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
258  int fieldNum = fieldIds_[fieldIndex];
259  const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
260 
261  // loop over basis functions
262  for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
263  int offset = elmtOffset[basis];
264  LO lid = LIDs[offset];
265  ScalarT value = (scatterFields_[fieldIndex])(worksetCellIndex,basis);
266  for(int d=0;d<value.size();d++)
267  dfdp_vectors_[d][lid] += value.fastAccessDx(d);
268  }
269  }
270  }
271 }
272 
273 // **********************************************************************
274 // Specialization: Jacobian
275 // **********************************************************************
276 
277 template<typename TRAITS,typename LO,typename GO,typename NodeT>
280  const Teuchos::ParameterList& p)
281  : globalIndexer_(indexer)
282  , globalDataKey_("Residual Scatter Container")
283 {
284  std::string scatterName = p.get<std::string>("Scatter Name");
285  scatterHolder_ =
286  Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
287 
288  // get names to be evaluated
289  const std::vector<std::string>& names =
290  *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
291 
292  // grab map from evaluated names to field names
293  fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
294 
296  p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional;
297 
298  // build the vector of fields that this is dependent on
299  scatterFields_.resize(names.size());
300  scratch_offsets_.resize(names.size());
301  for (std::size_t eq = 0; eq < names.size(); ++eq) {
302  scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
303 
304  // tell the field manager that we depend on this field
305  this->addDependentField(scatterFields_[eq]);
306  }
307 
308  // this is what this evaluator provides
309  this->addEvaluatedField(*scatterHolder_);
310 
311  if (p.isType<std::string>("Global Data Key"))
312  globalDataKey_ = p.get<std::string>("Global Data Key");
313 
314  this->setName(scatterName+" Scatter Residual (Jacobian)");
315 }
316 
317 // **********************************************************************
318 template<typename TRAITS,typename LO,typename GO,typename NodeT>
320 postRegistrationSetup(typename TRAITS::SetupData d,
321  PHX::FieldManager<TRAITS>& /* fm */)
322 {
323  fieldIds_.resize(scatterFields_.size());
324 
325  const Workset & workset_0 = (*d.worksets_)[0];
326  std::string blockId = this->wda(workset_0).block_id;
327 
328  // load required field numbers for fast use
329  for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
330  // get field ID from DOF manager
331  std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
332  fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
333 
334  int fieldNum = fieldIds_[fd];
335  const std::vector<int> & offsets = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
336  scratch_offsets_[fd] = Kokkos::View<int*,PHX::Device>("offsets",offsets.size());
337  for(std::size_t i=0;i<offsets.size();i++)
338  scratch_offsets_[fd](i) = offsets[i];
339  }
340 
341  scratch_lids_ = Kokkos::View<LO**,PHX::Device>("lids",scatterFields_[0].extent(0),
342  globalIndexer_->getElementBlockGIDCount(blockId));
343 }
344 
345 // **********************************************************************
346 template<typename TRAITS,typename LO,typename GO,typename NodeT>
348 preEvaluate(typename TRAITS::PreEvalData d)
349 {
351 
352  // extract linear object container
353  tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(d.gedc->getDataObject(globalDataKey_));
354 
355  if(tpetraContainer_==Teuchos::null) {
356  // extract linear object container
357  Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),true)->getGhostedLOC();
358  tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(loc);
359  }
360 }
361 
362 
363 // **********************************************************************
364 namespace panzer {
365 namespace {
366 
367 template <typename ScalarT,typename LO,typename GO,typename NodeT,typename LocalMatrixT>
368 class ScatterResidual_Jacobian_Functor {
369 public:
370  typedef typename PHX::Device execution_space;
371  typedef PHX::MDField<const ScalarT,Cell,NODE> FieldType;
372 
374  Kokkos::View<double**, Kokkos::LayoutLeft,PHX::Device> r_data;
375  LocalMatrixT jac; // Kokkos jacobian type
376 
377  Kokkos::View<const LO**,PHX::Device> lids; // local indices for unknowns
378  Kokkos::View<const int*,PHX::Device> offsets; // how to get a particular field
380 
381 
382  KOKKOS_INLINE_FUNCTION
383  void operator()(const unsigned int cell) const
384  {
385  LO cLIDs[256];
386  typename Sacado::ScalarType<ScalarT>::type vals[256];
387  int numIds = lids.extent(1);
388 
389  for(int i=0;i<numIds;i++)
390  cLIDs[i] = lids(cell,i);
391 
392  // loop over the basis functions (currently they are nodes)
393  for(std::size_t basis=0; basis < offsets.extent(0); basis++) {
394  typename FieldType::array_type::reference_type scatterField = field(cell,basis);
395  int offset = offsets(basis);
396  LO lid = lids(cell,offset);
397 
398  // Sum residual
399  if(fillResidual)
400  Kokkos::atomic_add(&r_data(lid,0), scatterField.val());
401 
402  // loop over the sensitivity indices: all DOFs on a cell
403  for(int sensIndex=0;sensIndex<numIds;++sensIndex)
404  vals[sensIndex] = scatterField.fastAccessDx(sensIndex);
405 
406  // Sum Jacobian
407  jac.sumIntoValues(lid, cLIDs,numIds, vals, true, true);
408  } // end basis
409  }
410 };
411 
412 template <typename ScalarT,typename LO,typename GO,typename NodeT>
413 class ScatterResidual_Residual_Functor {
414 public:
415  typedef typename PHX::Device execution_space;
416  typedef PHX::MDField<const ScalarT,Cell,NODE> FieldType;
417 
418  Kokkos::View<double**, Kokkos::LayoutLeft,PHX::Device> r_data;
419 
420  Kokkos::View<const LO**,PHX::Device> lids; // local indices for unknowns
421  Kokkos::View<const int*,PHX::Device> offsets; // how to get a particular field
423 
424 
425  KOKKOS_INLINE_FUNCTION
426  void operator()(const unsigned int cell) const
427  {
428 
429  // loop over the basis functions (currently they are nodes)
430  for(std::size_t basis=0; basis < offsets.extent(0); basis++) {
431  int offset = offsets(basis);
432  LO lid = lids(cell,offset);
433  Kokkos::atomic_add(&r_data(lid,0), field(cell,basis));
434 
435  } // end basis
436  }
437 };
438 
439 }
440 }
441 
442 // **********************************************************************
443 template<typename TRAITS,typename LO,typename GO,typename NodeT>
445 evaluateFields(typename TRAITS::EvalData workset)
446 {
448 
449  // for convenience pull out some objects from workset
450  std::string blockId = this->wda(workset).block_id;
451 
452  Teuchos::RCP<typename LOC::VectorType> r = tpetraContainer_->get_f();
453 
454  globalIndexer_->getElementLIDs(this->wda(workset).cell_local_ids_k,scratch_lids_);
455 
456  ScatterResidual_Residual_Functor<ScalarT,LO,GO,NodeT> functor;
457  functor.r_data = r->template getLocalView<PHX::Device>();
458  functor.lids = scratch_lids_;
459 
460  // for each field, do a parallel for loop
461  for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
462  functor.offsets = scratch_offsets_[fieldIndex];
463  functor.field = scatterFields_[fieldIndex];
464 
465  Kokkos::parallel_for(workset.num_cells,functor);
466  }
467 }
468 
469 // **********************************************************************
470 template<typename TRAITS,typename LO,typename GO,typename NodeT>
472 evaluateFields(typename TRAITS::EvalData workset)
473 {
475 
476  typedef typename LOC::CrsMatrixType::local_matrix_type LocalMatrixT;
477 
478  // for convenience pull out some objects from workset
479  std::string blockId = this->wda(workset).block_id;
480 
481  Teuchos::RCP<typename LOC::VectorType> r = tpetraContainer_->get_f();
482  Teuchos::RCP<typename LOC::CrsMatrixType> Jac = tpetraContainer_->get_A();
483 
484  globalIndexer_->getElementLIDs(this->wda(workset).cell_local_ids_k,scratch_lids_);
485 
486  ScatterResidual_Jacobian_Functor<ScalarT,LO,GO,NodeT,LocalMatrixT> functor;
487  functor.fillResidual = (r!=Teuchos::null);
488  if(functor.fillResidual)
489  functor.r_data = r->template getLocalView<PHX::Device>();
490  functor.jac = Jac->getLocalMatrix();
491  functor.lids = scratch_lids_;
492 
493  // for each field, do a parallel for loop
494  for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
495  functor.offsets = scratch_offsets_[fieldIndex];
496  functor.field = scatterFields_[fieldIndex];
497 
498  Kokkos::parallel_for(workset.num_cells,functor);
499  }
500 }
501 
502 // **********************************************************************
503 
504 #endif
T & get(const std::string &name, T def_value)
FieldType
The type of discretization to use for a field pattern.
Kokkos::View< const LO **, PHX::Device > lids
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Kokkos::View< double **, Kokkos::LayoutLeft, PHX::Device > r_data
Pushes residual values into the residual vector for a Newton-based solve.
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...
bool isType(const std::string &name) const
Kokkos::View< const int *, PHX::Device > offsets