Panzer  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Panzer_ScatterResidual_Epetra_impl.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Panzer: A partial differential equation assembly
4 // engine for strongly coupled complex multiphysics systems
5 //
6 // Copyright 2011 NTESS and the Panzer contributors.
7 // SPDX-License-Identifier: BSD-3-Clause
8 // *****************************************************************************
9 // @HEADER
10 
11 #ifndef PANZER_SCATTER_RESIDUAL_EPETRA_IMPL_HPP
12 #define PANZER_SCATTER_RESIDUAL_EPETRA_IMPL_HPP
13 
14 #include "Teuchos_RCP.hpp"
15 #include "Teuchos_Assert.hpp"
16 
17 #include "Phalanx_DataLayout.hpp"
18 
19 #include "Epetra_Map.h"
20 #include "Epetra_Vector.h"
21 #include "Epetra_CrsMatrix.h"
22 
23 #include "Panzer_GlobalIndexer.hpp"
24 #include "Panzer_PureBasis.hpp"
30 
31 #include "Phalanx_DataLayout_MDALayout.hpp"
32 
33 #include "Teuchos_FancyOStream.hpp"
34 
35 // **********************************************************************
36 // Specialization: Residual
37 // **********************************************************************
38 
39 template<typename TRAITS,typename LO,typename GO>
42  const Teuchos::RCP<const panzer::GlobalIndexer> & /* cIndexer */,
43  const Teuchos::ParameterList& p,
44  bool useDiscreteAdjoint)
45  : globalIndexer_(indexer)
46  , globalDataKey_("Residual Scatter Container")
47  , useDiscreteAdjoint_(useDiscreteAdjoint)
48 {
49  std::string scatterName = p.get<std::string>("Scatter Name");
50  scatterHolder_ =
52 
53  // get names to be evaluated
54  const std::vector<std::string>& names =
55  *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
56 
57  // grab map from evaluated names to field names
58  fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
59 
61  p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional;
62 
63  // build the vector of fields that this is dependent on
64  scatterFields_.resize(names.size());
65  for (std::size_t eq = 0; eq < names.size(); ++eq) {
66  scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
67 
68  // tell the field manager that we depend on this field
69  this->addDependentField(scatterFields_[eq]);
70  }
71 
72  // this is what this evaluator provides
73  this->addEvaluatedField(*scatterHolder_);
74 
75  if (p.isType<std::string>("Global Data Key"))
76  globalDataKey_ = p.get<std::string>("Global Data Key");
77 
78  this->setName(scatterName+" Scatter Residual");
79 }
80 
81 // **********************************************************************
82 template<typename TRAITS,typename LO,typename GO>
84 postRegistrationSetup(typename TRAITS::SetupData /* d */,
85  PHX::FieldManager<TRAITS>& /* fm */)
86 {
87  fieldIds_.resize(scatterFields_.size());
88  // load required field numbers for fast use
89  for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
90  // get field ID from DOF manager
91  std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
92  fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
93  }
94 }
95 
96 // **********************************************************************
97 template<typename TRAITS,typename LO,typename GO>
99 preEvaluate(typename TRAITS::PreEvalData d)
100 {
101  // extract linear object container
102  epetraContainer_ = Teuchos::rcp_dynamic_cast<EpetraLinearObjContainer>(d.gedc->getDataObject(globalDataKey_));
103 
104  if(epetraContainer_==Teuchos::null) {
105  // extract linear object container
106  Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),true)->getGhostedLOC();
107  epetraContainer_ = Teuchos::rcp_dynamic_cast<EpetraLinearObjContainer>(loc);
108  }
109 }
110 
111 // **********************************************************************
112 template<typename TRAITS,typename LO,typename GO>
114 evaluateFields(typename TRAITS::EvalData workset)
115 {
116  // for convenience pull out some objects from workset
117  std::string blockId = this->wda(workset).block_id;
118  const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
119 
120  Teuchos::RCP<Epetra_Vector> r = epetraContainer_->get_f();
121 
122  // NOTE: A reordering of these loops will likely improve performance
123  // The "getGIDFieldOffsets may be expensive. However the
124  // "getElementGIDs" can be cheaper. However the lookup for LIDs
125  // may be more expensive!
126 
127  // scatter operation for each cell in workset
128  auto LIDs = globalIndexer_->getLIDs();
129  auto LIDs_h = Kokkos::create_mirror_view(LIDs);
130  Kokkos::deep_copy(LIDs_h, LIDs);
131  // loop over each field to be scattered
132  for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
133  int fieldNum = fieldIds_[fieldIndex];
134  auto field = PHX::as_view(scatterFields_[fieldIndex]);
135  auto field_h = Kokkos::create_mirror_view(field);
136  Kokkos::deep_copy(field_h, field);
137 
138  const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
139  for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
140  std::size_t cellLocalId = localCellIds[worksetCellIndex];
141 
142  // loop over basis functions
143  for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
144  int offset = elmtOffset[basis];
145  int lid = LIDs_h(cellLocalId, offset);
146  (*r)[lid] += field_h(worksetCellIndex,basis);
147  }
148  }
149  }
150 }
151 
152 // **********************************************************************
153 // Specialization: Tangent
154 // **********************************************************************
155 
156 template<typename TRAITS,typename LO,typename GO>
159  const Teuchos::RCP<const panzer::GlobalIndexer> & /* cIndexer */,
160  const Teuchos::ParameterList& p,
161  bool /* useDiscreteAdjoint */)
162  : globalIndexer_(indexer)
163 {
164  std::string scatterName = p.get<std::string>("Scatter Name");
165  scatterHolder_ =
167 
168  // get names to be evaluated
169  const std::vector<std::string>& names =
170  *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
171 
172  // grab map from evaluated names to field names
173  fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
174 
176  p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional;
177 
178  // build the vector of fields that this is dependent on
179  scatterFields_.resize(names.size());
180  for (std::size_t eq = 0; eq < names.size(); ++eq) {
181  scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
182 
183  // tell the field manager that we depend on this field
184  this->addDependentField(scatterFields_[eq]);
185  }
186 
187  // this is what this evaluator provides
188  this->addEvaluatedField(*scatterHolder_);
189 
190  this->setName(scatterName+" Scatter Tangent");
191 }
192 
193 // **********************************************************************
194 template<typename TRAITS,typename LO,typename GO>
196 postRegistrationSetup(typename TRAITS::SetupData /* d */,
197  PHX::FieldManager<TRAITS>& /* fm */)
198 {
199  fieldIds_.resize(scatterFields_.size());
200  // load required field numbers for fast use
201  for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
202  // get field ID from DOF manager
203  std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
204  fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
205  }
206 }
207 
208 // **********************************************************************
209 template<typename TRAITS,typename LO,typename GO>
211 preEvaluate(typename TRAITS::PreEvalData d)
212 {
213  using Teuchos::RCP;
214  using Teuchos::rcp_dynamic_cast;
215 
216  // this is the list of parameters and their names that this scatter has to account for
217  std::vector<std::string> activeParameters =
218  rcp_dynamic_cast<ParameterList_GlobalEvaluationData>(d.gedc->getDataObject("PARAMETER_NAMES"))->getActiveParameters();
219 
220  dfdp_vectors_.clear();
221  for(std::size_t i=0;i<activeParameters.size();i++) {
222  RCP<Epetra_Vector> vec = rcp_dynamic_cast<EpetraLinearObjContainer>(d.gedc->getDataObject(activeParameters[i]),true)->get_f();
223  dfdp_vectors_.push_back(vec);
224  }
225 }
226 
227 // **********************************************************************
228 template<typename TRAITS,typename LO,typename GO>
230 evaluateFields(typename TRAITS::EvalData workset)
231 {
232  // for convenience pull out some objects from workset
233  std::string blockId = this->wda(workset).block_id;
234  const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
235 
236  // NOTE: A reordering of these loops will likely improve performance
237  // The "getGIDFieldOffsets may be expensive. However the
238  // "getElementGIDs" can be cheaper. However the lookup for LIDs
239  // may be more expensive!
240 
241  // loop over each field to be scattered
242  auto LIDs = globalIndexer_->getLIDs();
243  auto LIDs_h = Kokkos::create_mirror_view(LIDs);
244  Kokkos::deep_copy(LIDs_h, LIDs);
245  for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
246  int fieldNum = fieldIds_[fieldIndex];
247  const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
248  auto scatterField_h = Kokkos::create_mirror_view(scatterFields_[fieldIndex].get_static_view());
249  Kokkos::deep_copy(scatterField_h, scatterFields_[fieldIndex].get_static_view());
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  // loop over basis functions
255  for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
256  int offset = elmtOffset[basis];
257  int lid = LIDs_h(cellLocalId, offset);
258 
259  // scatter the sensitivity vectors
260  ScalarT value = scatterField_h(worksetCellIndex,basis);
261  for(int d=0;d<value.size();d++)
262  (*dfdp_vectors_[d])[lid] += value.fastAccessDx(d);
263  }
264  }
265  }
266 }
267 
268 // **********************************************************************
269 // Specialization: Jacobian
270 // **********************************************************************
271 
272 template<typename TRAITS,typename LO,typename GO>
276  const Teuchos::ParameterList& p,
277  bool useDiscreteAdjoint)
278  : globalIndexer_(indexer)
279  , colGlobalIndexer_(cIndexer)
280  , globalDataKey_("Residual Scatter Container")
281  , useDiscreteAdjoint_(useDiscreteAdjoint)
282 {
283  std::string scatterName = p.get<std::string>("Scatter Name");
284  scatterHolder_ =
286 
287  // get names to be evaluated
288  const std::vector<std::string>& names =
289  *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
290 
291  // grab map from evaluated names to field names
292  fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
293 
295  p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional;
296 
297  // build the vector of fields that this is dependent on
298  scatterFields_.resize(names.size());
299  for (std::size_t eq = 0; eq < names.size(); ++eq) {
300  scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
301 
302  // tell the field manager that we depend on this field
303  this->addDependentField(scatterFields_[eq]);
304  }
305 
306  // this is what this evaluator provides
307  this->addEvaluatedField(*scatterHolder_);
308 
309  if (p.isType<std::string>("Global Data Key"))
310  globalDataKey_ = p.get<std::string>("Global Data Key");
311  if (p.isType<bool>("Use Discrete Adjoint"))
312  useDiscreteAdjoint = p.get<bool>("Use Discrete Adjoint");
313 
314  // discrete adjoint does not work with non-square matrices
315  if(useDiscreteAdjoint)
316  { TEUCHOS_ASSERT(colGlobalIndexer_==globalIndexer_); }
317 
318  this->setName(scatterName+" Scatter Residual Epetra (Jacobian)");
319 }
320 
321 // **********************************************************************
322 template<typename TRAITS,typename LO,typename GO>
324 postRegistrationSetup(typename TRAITS::SetupData /* d */,
325  PHX::FieldManager<TRAITS>& /* fm */)
326 {
327  // globalIndexer_ = d.globalIndexer_;
328 
329  fieldIds_.resize(scatterFields_.size());
330  // load required field numbers for fast use
331  for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
332  // get field ID from DOF manager
333  std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
334  fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
335  }
336 }
337 
338 // **********************************************************************
339 template<typename TRAITS,typename LO,typename GO>
341 preEvaluate(typename TRAITS::PreEvalData d)
342 {
343  // extract linear object container
344  epetraContainer_ = Teuchos::rcp_dynamic_cast<EpetraLinearObjContainer>(d.gedc->getDataObject(globalDataKey_));
345 
346  if(epetraContainer_==Teuchos::null) {
347  // extract linear object container
348  Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),true)->getGhostedLOC();
349  epetraContainer_ = Teuchos::rcp_dynamic_cast<EpetraLinearObjContainer>(loc);
350  }
351 }
352 
353 // **********************************************************************
354 template<typename TRAITS,typename LO,typename GO>
356 evaluateFields(typename TRAITS::EvalData workset)
357 {
358  std::vector<double> jacRow;
359 
360  bool useColumnIndexer = colGlobalIndexer_!=Teuchos::null;
361 
362  // for convenience pull out some objects from workset
363  std::string blockId = this->wda(workset).block_id;
364  const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
365 
366  Teuchos::RCP<Epetra_Vector> r = epetraContainer_->get_f();
367  Teuchos::RCP<Epetra_CrsMatrix> Jac = epetraContainer_->get_A();
368 
370  colGlobalIndexer = useColumnIndexer ? colGlobalIndexer_ : globalIndexer_;
371 
372  // NOTE: A reordering of these loops will likely improve performance
373  // The "getGIDFieldOffsets" may be expensive. However the
374  // "getElementGIDs" can be cheaper. However the lookup for LIDs
375  // may be more expensive!
376 
377  // scatter operation for each cell in workset
378  auto LIDs = globalIndexer_->getLIDs();
379  auto colLIDs = colGlobalIndexer->getLIDs();
380  auto LIDs_h = Kokkos::create_mirror_view(LIDs);
381  auto colLIDs_h = Kokkos::create_mirror_view(colLIDs);
382  Kokkos::deep_copy(LIDs_h, LIDs);
383  Kokkos::deep_copy(colLIDs_h, colLIDs);
384  std::vector<typename decltype(scatterFields_[0].get_static_view())::host_mirror_type> scatterFields_h;
385  for ( std::size_t i=0; i< scatterFields_.size(); ++i) {
386  scatterFields_h.push_back(Kokkos::create_mirror_view(scatterFields_[i].get_static_view()));
387  Kokkos::deep_copy(scatterFields_h[i], scatterFields_[i].get_static_view());
388  }
389 
390  for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
391  std::size_t cellLocalId = localCellIds[worksetCellIndex];
392 
393  std::vector<int> cLIDs;
394  for (int i(0); i < static_cast<int>(colLIDs_h.extent(1)); ++i)
395  cLIDs.push_back(colLIDs_h(cellLocalId, i));
396  if (Teuchos::nonnull(workset.other)) {
397  const std::size_t other_cellLocalId = workset.other->cell_local_ids[worksetCellIndex];
398  for (int i(0); i < static_cast<int>(colLIDs_h.extent(1)); ++i)
399  cLIDs.push_back(colLIDs_h(other_cellLocalId, i));
400  }
401 
402  // loop over each field to be scattered
403  for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
404  int fieldNum = fieldIds_[fieldIndex];
405  const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
406 
407  // loop over the basis functions (currently they are nodes)
408  for(std::size_t rowBasisNum = 0; rowBasisNum < elmtOffset.size(); rowBasisNum++) {
409  const ScalarT scatterField = (scatterFields_h[fieldIndex])(worksetCellIndex,rowBasisNum);
410  int rowOffset = elmtOffset[rowBasisNum];
411  int row = LIDs_h(cellLocalId,rowOffset);
412 
413  // Sum residual
414  if(r!=Teuchos::null)
415  r->SumIntoMyValue(row,0,scatterField.val());
416 
417  // Sum Jacobian
418  if(useDiscreteAdjoint_) {
419  // loop over the sensitivity indices: all DOFs on a cell
420  jacRow.resize(scatterField.size());
421 
422  for(int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
423  jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex);
424 
425  for(std::size_t c=0;c<cLIDs.size();c++) {
426  int err = Jac->SumIntoMyValues(cLIDs[c], 1, &jacRow[c],&row);
428  }
429  }
430  else {
431  int err = Jac->SumIntoMyValues(
432  row,
433  std::min(cLIDs.size(), static_cast<size_t>(scatterField.size())),
434  scatterField.dx(),
435  panzer::ptrFromStlVector(cLIDs));
436 
438  }
439  } // end rowBasisNum
440  } // end fieldIndex
441  }
442 }
443 
444 // **********************************************************************
445 
446 #endif
T & get(const std::string &name, T def_value)
bool nonnull(const std::shared_ptr< T > &p)
const Kokkos::View< const panzer::LocalOrdinal **, Kokkos::LayoutRight, PHX::Device > getLIDs() const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
T * ptrFromStlVector(std::vector< T > &v)
Pushes residual values into the residual vector for a Newton-based solve.
int SumIntoMyValues(int MyRow, int NumEntries, const double *Values, const int *Indices)
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
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_ASSERT_EQUALITY(val1, val2)