Panzer  Version of the Day
Panzer_ScatterResidual_BlockedTpetra_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_BLOCKEDEPETRA_IMPL_HPP
44 #define PANZER_SCATTER_RESIDUAL_BLOCKEDEPETRA_IMPL_HPP
45 
46 #include "Teuchos_RCP.hpp"
47 #include "Teuchos_Assert.hpp"
48 
49 #include "Phalanx_DataLayout.hpp"
50 
53 #include "Panzer_PureBasis.hpp"
56 #include "Panzer_HashUtils.hpp"
57 
58 #include "Thyra_SpmdVectorBase.hpp"
59 #include "Thyra_ProductVectorBase.hpp"
60 #include "Thyra_BlockedLinearOpBase.hpp"
61 
62 #include "Phalanx_DataLayout_MDALayout.hpp"
63 
64 #include "Teuchos_FancyOStream.hpp"
65 
66 template <typename EvalT,typename TRAITS,typename LO,typename GO,typename NodeT>
69  const Teuchos::ParameterList& p)
70 {
71  std::string scatterName = p.get<std::string>("Scatter Name");
72  Teuchos::RCP<PHX::FieldTag> scatterHolder =
73  Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
74 
75  // get names to be evaluated
76  const std::vector<std::string>& names =
77  *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
78 
80  p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional;
81 
82  // build the vector of fields that this is dependent on
83  for (std::size_t eq = 0; eq < names.size(); ++eq) {
84  PHX::MDField<const ScalarT,Cell,NODE> field = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
85 
86  // tell the field manager that we depend on this field
87  this->addDependentField(field.fieldTag());
88  }
89 
90  // this is what this evaluator provides
91  this->addEvaluatedField(*scatterHolder);
92 
93  this->setName(scatterName+" Scatter Residual");
94 }
95 
96 // **********************************************************************
97 // Specialization: Residual
98 // **********************************************************************
99 
100 template <typename TRAITS,typename LO,typename GO,typename NodeT>
103  const Teuchos::ParameterList& p)
104  : globalIndexer_(indexer)
105  , globalDataKey_("Residual Scatter Container")
106 {
107  std::string scatterName = p.get<std::string>("Scatter Name");
108  scatterHolder_ =
109  Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
110 
111  // get names to be evaluated
112  const std::vector<std::string>& names =
113  *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
114 
115  // grab map from evaluated names to field names
116  fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
117 
119  p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional;
120 
121  // build the vector of fields that this is dependent on
122  scatterFields_.resize(names.size());
123  for (std::size_t eq = 0; eq < names.size(); ++eq) {
124  scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
125 
126  // tell the field manager that we depend on this field
127  this->addDependentField(scatterFields_[eq]);
128  }
129 
130  // this is what this evaluator provides
131  this->addEvaluatedField(*scatterHolder_);
132 
133  if (p.isType<std::string>("Global Data Key"))
134  globalDataKey_ = p.get<std::string>("Global Data Key");
135 
136  this->setName(scatterName+" Scatter Residual");
137 }
138 
139 // **********************************************************************
140 template <typename TRAITS,typename LO,typename GO,typename NodeT>
142 postRegistrationSetup(typename TRAITS::SetupData d,
144 {
145  fieldIds_.resize(scatterFields_.size());
146 
147  // load required field numbers for fast use
148  for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
149  // get field ID from DOF manager
150  std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
151  fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
152 
153  // fill field data object
154  this->utils.setFieldData(scatterFields_[fd],fm);
155  }
156 }
157 
158 // **********************************************************************
159 template <typename TRAITS,typename LO,typename GO,typename NodeT>
161 preEvaluate(typename TRAITS::PreEvalData d)
162 {
163  // extract linear object container
164  blockedContainer_ = Teuchos::rcp_dynamic_cast<const ContainerType>(d.gedc.getDataObject(globalDataKey_),true);
165 }
166 
167 // **********************************************************************
168 template <typename TRAITS,typename LO,typename GO,typename NodeT>
170 evaluateFields(typename TRAITS::EvalData workset)
171 {
172  using Teuchos::RCP;
173  using Teuchos::ArrayRCP;
174  using Teuchos::ptrFromRef;
175  using Teuchos::rcp_dynamic_cast;
176 
177  using Thyra::VectorBase;
178  using Thyra::SpmdVectorBase;
180 
181  std::vector<std::pair<int,GO> > GIDs;
182  std::vector<LO> LIDs;
183 
184  // for convenience pull out some objects from workset
185  std::string blockId = this->wda(workset).block_id;
186  const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
187 
188  RCP<const ContainerType> blockedContainer = blockedContainer_;
189 
190  RCP<ProductVectorBase<double> > r = rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer->get_f(),true);
191 
192  // NOTE: A reordering of these loops will likely improve performance
193  // The "getGIDFieldOffsets may be expensive. However the
194  // "getElementGIDs" can be cheaper. However the lookup for LIDs
195  // may be more expensive!
196 
197  // scatter operation for each cell in workset
198  for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
199  std::size_t cellLocalId = localCellIds[worksetCellIndex];
200 
201  globalIndexer_->getElementGIDs(cellLocalId,GIDs,blockId);
202 
203  // caculate the local IDs for this element
204  LIDs.resize(GIDs.size());
205  for(std::size_t i=0;i<GIDs.size();i++) {
206  // used for doing local ID lookups
207  RCP<const MapType> r_map = blockedContainer->getMapForBlock(GIDs[i].first);
208 
209  LIDs[i] = r_map->getLocalElement(GIDs[i].second);
210  }
211 
212  // loop over each field to be scattered
214  for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
215  int fieldNum = fieldIds_[fieldIndex];
216  int indexerId = globalIndexer_->getFieldBlock(fieldNum);
217 
218  // grab local data for inputing
219  RCP<SpmdVectorBase<double> > block_r = rcp_dynamic_cast<SpmdVectorBase<double> >(r->getNonconstVectorBlock(indexerId));
220  block_r->getNonconstLocalData(ptrFromRef(local_r));
221 
222  const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
223 
224  // loop over basis functions
225  for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
226  int offset = elmtOffset[basis];
227  int lid = LIDs[offset];
228  local_r[lid] += (scatterFields_[fieldIndex])(worksetCellIndex,basis);
229  }
230  }
231  }
232 }
233 
234 // **********************************************************************
235 // Specialization: Jacobian
236 // **********************************************************************
237 
238 template <typename TRAITS,typename LO,typename GO,typename NodeT>
241  const Teuchos::ParameterList& p)
242  : globalIndexer_(indexer)
243  , globalDataKey_("Residual Scatter Container")
244 {
245  std::string scatterName = p.get<std::string>("Scatter Name");
246  scatterHolder_ =
247  Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
248 
249  // get names to be evaluated
250  const std::vector<std::string>& names =
251  *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
252 
253  // grab map from evaluated names to field names
254  fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
255 
257  p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional;
258 
259  // build the vector of fields that this is dependent on
260  scatterFields_.resize(names.size());
261  for (std::size_t eq = 0; eq < names.size(); ++eq) {
262  scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
263 
264  // tell the field manager that we depend on this field
265  this->addDependentField(scatterFields_[eq]);
266  }
267 
268  // this is what this evaluator provides
269  this->addEvaluatedField(*scatterHolder_);
270 
271  if (p.isType<std::string>("Global Data Key"))
272  globalDataKey_ = p.get<std::string>("Global Data Key");
273 
274  this->setName(scatterName+" Scatter Residual (Jacobian)");
275 }
276 
277 // **********************************************************************
278 template <typename TRAITS,typename LO,typename GO,typename NodeT>
280 postRegistrationSetup(typename TRAITS::SetupData d,
282 {
283  fieldIds_.resize(scatterFields_.size());
284 
285  // load required field numbers for fast use
286  for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
287  // get field ID from DOF manager
288  std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
289  fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
290 
291  // fill field data object
292  this->utils.setFieldData(scatterFields_[fd],fm);
293  }
294 }
295 
296 // **********************************************************************
297 template <typename TRAITS,typename LO,typename GO,typename NodeT>
299 preEvaluate(typename TRAITS::PreEvalData d)
300 {
301  using Teuchos::RCP;
302  using Teuchos::rcp_dynamic_cast;
303 
304  // extract linear object container
305  blockedContainer_ = rcp_dynamic_cast<const ContainerType>(d.gedc.getDataObject(globalDataKey_));
306 
307  if(blockedContainer_==Teuchos::null) {
308  RCP<const LOCPair_GlobalEvaluationData> gdata = rcp_dynamic_cast<const LOCPair_GlobalEvaluationData>(d.gedc.getDataObject(globalDataKey_),true);
309  blockedContainer_ = rcp_dynamic_cast<const ContainerType>(gdata->getGhostedLOC());
310  }
311 }
312 
313 // **********************************************************************
314 template <typename TRAITS,typename LO,typename GO,typename NodeT>
316 evaluateFields(typename TRAITS::EvalData workset)
317 {
318  using Teuchos::RCP;
319  using Teuchos::ArrayRCP;
320  using Teuchos::ptrFromRef;
321  using Teuchos::rcp_dynamic_cast;
322 
323  using Thyra::VectorBase;
324  using Thyra::SpmdVectorBase;
327 
328  std::vector<std::pair<int,GO> > GIDs;
329  std::vector<LO> LIDs;
330  std::vector<double> jacRow;
331 
332  // for convenience pull out some objects from workset
333  std::string blockId = this->wda(workset).block_id;
334  const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
335 
336  RCP<const ContainerType> blockedContainer = blockedContainer_;
337 
338  RCP<ProductVectorBase<double> > r = rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer->get_f());
339  Teuchos::RCP<BlockedLinearOpBase<double> > Jac = rcp_dynamic_cast<BlockedLinearOpBase<double> >(blockedContainer->get_A());
340 
341  int numFieldBlocks = globalIndexer_->getNumFieldBlocks();
342  std::vector<int> blockOffsets(numFieldBlocks+1); // number of fields, plus a sentinnel
343  for(int blk=0;blk<numFieldBlocks;blk++) {
344  int blockOffset = globalIndexer_->getBlockGIDOffset(blockId,blk);
345  blockOffsets[blk] = blockOffset;
346  }
347 
348  std::unordered_map<std::pair<int,int>,Teuchos::RCP<CrsMatrixType>,panzer::pair_hash> jacTpetraBlocks;
349 
350  // NOTE: A reordering of these loops will likely improve performance
351  // The "getGIDFieldOffsets" may be expensive. However the
352  // "getElementGIDs" can be cheaper. However the lookup for LIDs
353  // may be more expensive!
354 
355  // scatter operation for each cell in workset
356  for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
357  std::size_t cellLocalId = localCellIds[worksetCellIndex];
358 
359  globalIndexer_->getElementGIDs(cellLocalId,GIDs,blockId);
360 
361  // caculate the local IDs for this element
362  LIDs.resize(GIDs.size());
363  for(std::size_t i=0;i<GIDs.size();i++) {
364  // used for doing local ID lookups
365  RCP<const MapType> r_map = blockedContainer->getMapForBlock(GIDs[i].first);
366 
367  LIDs[i] = r_map->getLocalElement(GIDs[i].second);
368  }
369 
370  // loop over each field to be scattered
372  for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
373  int fieldNum = fieldIds_[fieldIndex];
374  int blockRowIndex = globalIndexer_->getFieldBlock(fieldNum);
375 
376  // grab local data for inputing
377  if(r!=Teuchos::null) {
378  RCP<SpmdVectorBase<double> > block_r = rcp_dynamic_cast<SpmdVectorBase<double> >(r->getNonconstVectorBlock(blockRowIndex));
379  block_r->getNonconstLocalData(ptrFromRef(local_r));
380  }
381 
382  const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
383 
384  // loop over the basis functions (currently they are nodes)
385  for(std::size_t rowBasisNum = 0; rowBasisNum < elmtOffset.size(); rowBasisNum++) {
386  const ScalarT scatterField = (scatterFields_[fieldIndex])(worksetCellIndex,rowBasisNum);
387  int rowOffset = elmtOffset[rowBasisNum];
388  int r_lid = LIDs[rowOffset];
389 
390  // Sum residual
391  if(local_r!=Teuchos::null)
392  local_r[r_lid] += (scatterField.val());
393 
394  blockOffsets[numFieldBlocks] = scatterField.size(); // add the sentinel
395 
396  // loop over the sensitivity indices: all DOFs on a cell
397  jacRow.resize(scatterField.size());
398 
399  for(int sensIndex=0;sensIndex<scatterField.size();++sensIndex) {
400  jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex);
401  }
402 
403  for(int blockColIndex=0;blockColIndex<numFieldBlocks;blockColIndex++) {
404  int start = blockOffsets[blockColIndex];
405  int end = blockOffsets[blockColIndex+1];
406 
407  if(end-start<=0)
408  continue;
409 
410  // check hash table for jacobian sub block
411  std::pair<int,int> blockIndex = std::make_pair(blockRowIndex,blockColIndex);
412  Teuchos::RCP<CrsMatrixType> subJac = jacTpetraBlocks[blockIndex];
413 
414  // if you didn't find one before, add it to the hash table
415  if(subJac==Teuchos::null) {
416  Teuchos::RCP<Thyra::LinearOpBase<double> > tOp = Jac->getNonconstBlock(blockIndex.first,blockIndex.second);
417 
418  // block operator is null, don't do anything (it is excluded)
419  if(Teuchos::is_null(tOp))
420  continue;
421 
422  Teuchos::RCP<OperatorType> tpetra_Op = rcp_dynamic_cast<ThyraLinearOp>(tOp)->getTpetraOperator();
423  subJac = rcp_dynamic_cast<CrsMatrixType>(tpetra_Op,true);
424  jacTpetraBlocks[blockIndex] = subJac;
425  }
426 
427  // Sum Jacobian
428  subJac->sumIntoLocalValues(r_lid, Teuchos::arrayViewFromVector(LIDs).view(start,end-start),
429  Teuchos::arrayViewFromVector(jacRow).view(start,end-start));
430  }
431  } // end rowBasisNum
432  } // end fieldIndex
433  }
434 }
435 
436 // **********************************************************************
437 
438 #endif
Pushes residual values into the residual vector for a Newton-based solve.
std::vector< PHX::MDField< const ScalarT, panzer::Cell, panzer::Point > > scatterFields_
size_type size() const
T & get(const std::string &name, T def_value)
bool is_null(const std::shared_ptr< T > &p)
ScatterResidual_BlockedTpetra(const Teuchos::RCP< const BlockedDOFManager< LO, GO > > &indexer)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
bool isType(const std::string &name) const
PHX::MDField< const ScalarT, Cell, IP > field
Teuchos::RCP< LinearObjContainer > getGhostedLOC() const
Teuchos::RCP< const panzer::PureBasis > basis
Interpolates basis DOF values to IP DOF values.
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &vm)