Panzer  Version of the Day
Panzer_ScatterResidual_Epetra_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_EPETRA_IMPL_HPP
44 #define PANZER_SCATTER_RESIDUAL_EPETRA_IMPL_HPP
45 
46 #include "Teuchos_RCP.hpp"
47 #include "Teuchos_Assert.hpp"
48 
49 #include "Phalanx_DataLayout.hpp"
50 
51 #include "Epetra_Map.h"
52 #include "Epetra_Vector.h"
53 #include "Epetra_CrsMatrix.h"
54 
56 #include "Panzer_PureBasis.hpp"
61 
62 #include "Phalanx_DataLayout_MDALayout.hpp"
63 
64 #include "Teuchos_FancyOStream.hpp"
65 
66 // **********************************************************************
67 // Specialization: Residual
68 // **********************************************************************
69 
70 template<typename TRAITS,typename LO,typename GO>
73  const Teuchos::RCP<const panzer::UniqueGlobalIndexer<LO,GO> > & cIndexer,
74  const Teuchos::ParameterList& p,
75  bool useDiscreteAdjoint)
76  : globalIndexer_(indexer)
77  , globalDataKey_("Residual Scatter Container")
78  , useDiscreteAdjoint_(useDiscreteAdjoint)
79 {
80  std::string scatterName = p.get<std::string>("Scatter Name");
81  scatterHolder_ =
82  Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
83 
84  // get names to be evaluated
85  const std::vector<std::string>& names =
86  *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
87 
88  // grab map from evaluated names to field names
89  fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
90 
92  p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional;
93 
94  // build the vector of fields that this is dependent on
95  scatterFields_.resize(names.size());
96  for (std::size_t eq = 0; eq < names.size(); ++eq) {
97  scatterFields_[eq] = PHX::MDField<ScalarT,Cell,NODE>(names[eq],dl);
98 
99  // tell the field manager that we depend on this field
100  this->addDependentField(scatterFields_[eq]);
101  }
102 
103  // this is what this evaluator provides
104  this->addEvaluatedField(*scatterHolder_);
105 
106  if (p.isType<std::string>("Global Data Key"))
107  globalDataKey_ = p.get<std::string>("Global Data Key");
108 
109  this->setName(scatterName+" Scatter Residual");
110 }
111 
112 // **********************************************************************
113 template<typename TRAITS,typename LO,typename GO>
115 postRegistrationSetup(typename TRAITS::SetupData d,
117 {
118  fieldIds_.resize(scatterFields_.size());
119  // load required field numbers for fast use
120  for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
121  // get field ID from DOF manager
122  std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
123  fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
124 
125  // fill field data object
126  this->utils.setFieldData(scatterFields_[fd],fm);
127  }
128 }
129 
130 // **********************************************************************
131 template<typename TRAITS,typename LO,typename GO>
133 preEvaluate(typename TRAITS::PreEvalData d)
134 {
135  // extract linear object container
136  epetraContainer_ = Teuchos::rcp_dynamic_cast<EpetraLinearObjContainer>(d.gedc.getDataObject(globalDataKey_));
137 
138  if(epetraContainer_==Teuchos::null) {
139  // extract linear object container
140  Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<LOCPair_GlobalEvaluationData>(d.gedc.getDataObject(globalDataKey_),true)->getGhostedLOC();
141  epetraContainer_ = Teuchos::rcp_dynamic_cast<EpetraLinearObjContainer>(loc);
142  }
143 }
144 
145 // **********************************************************************
146 template<typename TRAITS,typename LO,typename GO>
148 evaluateFields(typename TRAITS::EvalData workset)
149 {
150  std::vector<int> 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<Epetra_Vector> r = epetraContainer_->get_f();
157 
158  // NOTE: A reordering of these loops will likely improve performance
159  // The "getGIDFieldOffsets may be expensive. However the
160  // "getElementGIDs" can be cheaper. However the lookup for LIDs
161  // may be more expensive!
162 
163  // scatter operation for each cell in workset
164  for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
165  std::size_t cellLocalId = localCellIds[worksetCellIndex];
166 
167  LIDs = globalIndexer_->getElementLIDs(cellLocalId);
168 
169  // loop over each field to be scattered
170  for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
171  int fieldNum = fieldIds_[fieldIndex];
172  const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
173 
174  // loop over basis functions
175  for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
176  int offset = elmtOffset[basis];
177  int lid = LIDs[offset];
178  (*r)[lid] += (scatterFields_[fieldIndex])(worksetCellIndex,basis);
179  }
180  }
181  }
182 }
183 
184 // **********************************************************************
185 // Specialization: Tangent
186 // **********************************************************************
187 
188 template<typename TRAITS,typename LO,typename GO>
191  const Teuchos::RCP<const panzer::UniqueGlobalIndexer<LO,GO> > & cIndexer,
192  const Teuchos::ParameterList& p,
193  bool useDiscreteAdjoint)
194  : globalIndexer_(indexer)
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 
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<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  this->setName(scatterName+" Scatter Tangent");
223 }
224 
225 // **********************************************************************
226 template<typename TRAITS,typename LO,typename GO>
228 postRegistrationSetup(typename TRAITS::SetupData d,
230 {
231  fieldIds_.resize(scatterFields_.size());
232  // load required field numbers for fast use
233  for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
234  // get field ID from DOF manager
235  std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
236  fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
237 
238  // fill field data object
239  this->utils.setFieldData(scatterFields_[fd],fm);
240  }
241 }
242 
243 // **********************************************************************
244 template<typename TRAITS,typename LO,typename GO>
246 preEvaluate(typename TRAITS::PreEvalData d)
247 {
248  using Teuchos::RCP;
249  using Teuchos::rcp_dynamic_cast;
250 
251  // this is the list of parameters and their names that this scatter has to account for
252  std::vector<std::string> activeParameters =
253  rcp_dynamic_cast<ParameterList_GlobalEvaluationData>(d.gedc.getDataObject("PARAMETER_NAMES"))->getActiveParameters();
254 
255  dfdp_vectors_.clear();
256  for(std::size_t i=0;i<activeParameters.size();i++) {
257  RCP<Epetra_Vector> vec = rcp_dynamic_cast<EpetraLinearObjContainer>(d.gedc.getDataObject(activeParameters[i]),true)->get_f();
258  dfdp_vectors_.push_back(vec);
259  }
260 }
261 
262 // **********************************************************************
263 template<typename TRAITS,typename LO,typename GO>
265 evaluateFields(typename TRAITS::EvalData workset)
266 {
267  std::vector<int> LIDs;
268 
269  // for convenience pull out some objects from workset
270  std::string blockId = this->wda(workset).block_id;
271  const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
272 
273  // NOTE: A reordering of these loops will likely improve performance
274  // The "getGIDFieldOffsets may be expensive. However the
275  // "getElementGIDs" can be cheaper. However the lookup for LIDs
276  // may be more expensive!
277 
278  // scatter operation for each cell in workset
279  for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
280  std::size_t cellLocalId = localCellIds[worksetCellIndex];
281 
282  LIDs = globalIndexer_->getElementLIDs(cellLocalId);
283 
284  // loop over each field to be scattered
285  for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
286  int fieldNum = fieldIds_[fieldIndex];
287  const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
288 
289  // loop over basis functions
290  for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
291  int offset = elmtOffset[basis];
292  int lid = LIDs[offset];
293 
294  // scatter the sensitivity vectors
295  ScalarT value = (scatterFields_[fieldIndex])(worksetCellIndex,basis);
296  for(int d=0;d<value.size();d++)
297  (*dfdp_vectors_[d])[lid] += value.fastAccessDx(d);
298  }
299  }
300  }
301 }
302 
303 // **********************************************************************
304 // Specialization: Jacobian
305 // **********************************************************************
306 
307 template<typename TRAITS,typename LO,typename GO>
310  const Teuchos::RCP<const panzer::UniqueGlobalIndexer<LO,GO> > & cIndexer,
311  const Teuchos::ParameterList& p,
312  bool useDiscreteAdjoint)
313  : globalIndexer_(indexer)
314  , colGlobalIndexer_(cIndexer)
315  , globalDataKey_("Residual Scatter Container")
316  , useDiscreteAdjoint_(useDiscreteAdjoint)
317 {
318  std::string scatterName = p.get<std::string>("Scatter Name");
319  scatterHolder_ =
320  Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
321 
322  // get names to be evaluated
323  const std::vector<std::string>& names =
324  *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
325 
326  // grab map from evaluated names to field names
327  fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
328 
330  p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional;
331 
332  // build the vector of fields that this is dependent on
333  scatterFields_.resize(names.size());
334  for (std::size_t eq = 0; eq < names.size(); ++eq) {
335  scatterFields_[eq] = PHX::MDField<ScalarT,Cell,NODE>(names[eq],dl);
336 
337  // tell the field manager that we depend on this field
338  this->addDependentField(scatterFields_[eq]);
339  }
340 
341  // this is what this evaluator provides
342  this->addEvaluatedField(*scatterHolder_);
343 
344  if (p.isType<std::string>("Global Data Key"))
345  globalDataKey_ = p.get<std::string>("Global Data Key");
346  if (p.isType<bool>("Use Discrete Adjoint"))
347  useDiscreteAdjoint = p.get<bool>("Use Discrete Adjoint");
348 
349  // discrete adjoint does not work with non-square matrices
350  if(useDiscreteAdjoint)
351  { TEUCHOS_ASSERT(colGlobalIndexer_==globalIndexer_); }
352 
353  this->setName(scatterName+" Scatter Residual Epetra (Jacobian)");
354 }
355 
356 // **********************************************************************
357 template<typename TRAITS,typename LO,typename GO>
359 postRegistrationSetup(typename TRAITS::SetupData d,
361 {
362  // globalIndexer_ = d.globalIndexer_;
363 
364  fieldIds_.resize(scatterFields_.size());
365  // load required field numbers for fast use
366  for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
367  // get field ID from DOF manager
368  std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
369  fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
370 
371  // fill field data object
372  this->utils.setFieldData(scatterFields_[fd],fm);
373  }
374 }
375 
376 // **********************************************************************
377 template<typename TRAITS,typename LO,typename GO>
379 preEvaluate(typename TRAITS::PreEvalData d)
380 {
381  // extract linear object container
382  epetraContainer_ = Teuchos::rcp_dynamic_cast<EpetraLinearObjContainer>(d.gedc.getDataObject(globalDataKey_));
383 
384  if(epetraContainer_==Teuchos::null) {
385  // extract linear object container
386  Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<LOCPair_GlobalEvaluationData>(d.gedc.getDataObject(globalDataKey_),true)->getGhostedLOC();
387  epetraContainer_ = Teuchos::rcp_dynamic_cast<EpetraLinearObjContainer>(loc);
388  }
389 }
390 
391 // **********************************************************************
392 template<typename TRAITS,typename LO,typename GO>
394 evaluateFields(typename TRAITS::EvalData workset)
395 {
396  std::vector<int> cLIDs, rLIDs;
397  std::vector<double> jacRow;
398 
399  bool useColumnIndexer = colGlobalIndexer_!=Teuchos::null;
400 
401  // for convenience pull out some objects from workset
402  std::string blockId = this->wda(workset).block_id;
403  const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
404 
405  Teuchos::RCP<Epetra_Vector> r = epetraContainer_->get_f();
406  Teuchos::RCP<Epetra_CrsMatrix> Jac = epetraContainer_->get_A();
407 
409  colGlobalIndexer = useColumnIndexer ? colGlobalIndexer_ : globalIndexer_;
410 
411  // NOTE: A reordering of these loops will likely improve performance
412  // The "getGIDFieldOffsets" may be expensive. However the
413  // "getElementGIDs" can be cheaper. However the lookup for LIDs
414  // may be more expensive!
415 
416  // scatter operation for each cell in workset
417  for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
418  std::size_t cellLocalId = localCellIds[worksetCellIndex];
419 
420  rLIDs = globalIndexer_->getElementLIDs(cellLocalId);
421  cLIDs = colGlobalIndexer->getElementLIDs(cellLocalId);
422  if (Teuchos::nonnull(workset.other)) {
423  const std::size_t other_cellLocalId = workset.other->cell_local_ids[worksetCellIndex];
424  const std::vector<int> other_cLIDs = colGlobalIndexer->getElementLIDs(other_cellLocalId);
425  cLIDs.insert(cLIDs.end(), other_cLIDs.begin(), other_cLIDs.end());
426  }
427 
428  // loop over each field to be scattered
429  for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
430  int fieldNum = fieldIds_[fieldIndex];
431  const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
432 
433  // loop over the basis functions (currently they are nodes)
434  for(std::size_t rowBasisNum = 0; rowBasisNum < elmtOffset.size(); rowBasisNum++) {
435  const ScalarT scatterField = (scatterFields_[fieldIndex])(worksetCellIndex,rowBasisNum);
436  int rowOffset = elmtOffset[rowBasisNum];
437  int row = rLIDs[rowOffset];
438 
439  // Sum residual
440  if(r!=Teuchos::null)
441  r->SumIntoMyValue(row,0,scatterField.val());
442 
443  // Sum Jacobian
444  if(useDiscreteAdjoint_) {
445  // loop over the sensitivity indices: all DOFs on a cell
446  jacRow.resize(scatterField.size());
447 
448  for(int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
449  jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex);
450 
451  for(std::size_t c=0;c<cLIDs.size();c++) {
452  int err = Jac->SumIntoMyValues(cLIDs[c], 1, &jacRow[c],&row);
454  }
455  }
456  else {
457  int err = Jac->SumIntoMyValues(
458  row,
459  std::min(cLIDs.size(), static_cast<size_t>(scatterField.size())),
460  scatterField.dx(),
461  panzer::ptrFromStlVector(cLIDs));
463  }
464  } // end rowBasisNum
465  } // end fieldIndex
466  }
467 }
468 
469 // **********************************************************************
470 
471 #endif
std::vector< PHX::MDField< const ScalarT, panzer::Cell, panzer::Point > > scatterFields_
T & get(const std::string &name, T def_value)
bool nonnull(const std::shared_ptr< T > &p)
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 SumIntoMyValue(int MyRow, int VectorIndex, double ScalarValue)
bool isType(const std::string &name) const
int SumIntoMyValues(int MyRow, int NumEntries, const double *Values, const int *Indices)
Teuchos::RCP< const panzer::PureBasis > basis
Interpolates basis DOF values to IP DOF values.
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_ASSERT_EQUALITY(val1, val2)