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