Panzer  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Panzer_ScatterDirichletResidual_BlockedEpetra_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_DIRICHLET_RESIDUAL_BLOCEDEPETRA_IMPL_HPP
12 #define PANZER_SCATTER_DIRICHLET_RESIDUAL_BLOCEDEPETRA_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_GlobalIndexer_Utilities.hpp"
25 #include "Panzer_BlockedDOFManager.hpp"
26 #include "Panzer_PureBasis.hpp"
27 #include "Panzer_BlockedEpetraLinearObjContainer.hpp"
28 #include "Panzer_GlobalEvaluationDataContainer.hpp"
29 
30 #include "Phalanx_DataLayout_MDALayout.hpp"
31 
32 #include "Thyra_SpmdVectorBase.hpp"
33 #include "Thyra_ProductVectorBase.hpp"
34 #include "Thyra_DefaultProductVector.hpp"
35 #include "Thyra_BlockedLinearOpBase.hpp"
36 #include "Thyra_get_Epetra_Operator.hpp"
37 
38 #include "Teuchos_FancyOStream.hpp"
39 
40 #include <unordered_map>
41 
42 // **********************************************************************
43 // Specialization: Residual
44 // **********************************************************************
45 
46 
47 template<typename TRAITS,typename LO,typename GO>
48 panzer::ScatterDirichletResidual_BlockedEpetra<panzer::Traits::Residual, TRAITS,LO,GO>::
49 ScatterDirichletResidual_BlockedEpetra(const std::vector<Teuchos::RCP<const GlobalIndexer> > & rIndexers,
50  const std::vector<Teuchos::RCP<const GlobalIndexer> > & cIndexers,
51  const Teuchos::ParameterList& p,
52  bool /* useDiscreteAdjoint */)
53  : rowIndexers_(rIndexers)
54  , colIndexers_(cIndexers)
55  , globalDataKey_("Residual Scatter Container")
56 {
57  std::string scatterName = p.get<std::string>("Scatter Name");
58  scatterHolder_ =
59  Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
60 
61  // get names to be evaluated
62  const std::vector<std::string>& names =
63  *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
64 
65  // grab map from evaluated names to field names
66  fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
67 
68  // determine if we are scattering an initial condition
69  scatterIC_ = p.isParameter("Scatter Initial Condition") ? p.get<bool>("Scatter Initial Condition") : false;
70 
71  Teuchos::RCP<PHX::DataLayout> dl = (!scatterIC_) ?
72  p.get< Teuchos::RCP<panzer::PureBasis> >("Basis")->functional :
73  p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional;
74  if (!scatterIC_) {
75  side_subcell_dim_ = p.get<int>("Side Subcell Dimension");
76  local_side_id_ = p.get<int>("Local Side ID");
77  }
78 
79  // build the vector of fields that this is dependent on
80  scatterFields_.resize(names.size());
81  for (std::size_t eq = 0; eq < names.size(); ++eq) {
82  scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
83 
84  // tell the field manager that we depend on this field
85  this->addDependentField(scatterFields_[eq]);
86  }
87 
88  checkApplyBC_ = p.isParameter("Check Apply BC") ? p.get<bool>("Check Apply BC") : false;
89  if (checkApplyBC_) {
90  applyBC_.resize(names.size());
91  for (std::size_t eq = 0; eq < names.size(); ++eq) {
92  applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string("APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
93  this->addDependentField(applyBC_[eq]);
94  }
95  }
96 
97  // this is what this evaluator provides
98  this->addEvaluatedField(*scatterHolder_);
99 
100  if (p.isType<std::string>("Global Data Key"))
101  globalDataKey_ = p.get<std::string>("Global Data Key");
102 
103  this->setName(scatterName+" Scatter Residual");
104 }
105 
106 // **********************************************************************
107 template<typename TRAITS,typename LO,typename GO>
108 void panzer::ScatterDirichletResidual_BlockedEpetra<panzer::Traits::Residual, TRAITS,LO,GO>::
109 postRegistrationSetup(typename TRAITS::SetupData /* d */,
110  PHX::FieldManager<TRAITS>& /* fm */)
111 {
112  indexerIds_.resize(scatterFields_.size());
113  subFieldIds_.resize(scatterFields_.size());
114 
115  // load required field numbers for fast use
116  for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
117  // get field ID from DOF manager
118  std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
119 
120  indexerIds_[fd] = getFieldBlock(fieldName,rowIndexers_);
121  subFieldIds_[fd] = rowIndexers_[indexerIds_[fd]]->getFieldNum(fieldName);
122  }
123 
124  // get the number of nodes (Should be renamed basis)
125  num_nodes = scatterFields_[0].extent(1);
126 }
127 
128 // **********************************************************************
129 template<typename TRAITS,typename LO,typename GO>
130 void panzer::ScatterDirichletResidual_BlockedEpetra<panzer::Traits::Residual, TRAITS,LO,GO>::
131 preEvaluate(typename TRAITS::PreEvalData d)
132 {
133  typedef BlockedEpetraLinearObjContainer BLOC;
134  typedef BlockedEpetraLinearObjContainer ELOC;
135 
136  using Teuchos::rcp_dynamic_cast;
137  using Thyra::ProductVectorBase;
138 
139  // extract dirichlet counter from container
140  Teuchos::RCP<BLOC> blockContainer
141  = Teuchos::rcp_dynamic_cast<BLOC>(d.gedc->getDataObject("Dirichlet Counter"),true);
142 
143  dirichletCounter_ = Teuchos::rcp_dynamic_cast<Thyra::ProductVectorBase<double> >(blockContainer->get_f(),true);
144  TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
145 
146  // extract linear object container
147  Teuchos::RCP<const BLOC> blockedContainer = Teuchos::rcp_dynamic_cast<const BLOC>(d.gedc->getDataObject(globalDataKey_));
148  Teuchos::RCP<const ELOC> epetraContainer = Teuchos::rcp_dynamic_cast<const ELOC>(d.gedc->getDataObject(globalDataKey_));
149 
150  // if its blocked do this
151  if(blockedContainer!=Teuchos::null)
152  r_ = (!scatterIC_) ?
153  rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer->get_f(),true) :
154  rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer->get_x(),true);
155  else if(epetraContainer!=Teuchos::null) // if its straight up epetra do this
156  r_ = (!scatterIC_) ?
157  Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_f_th()) :
158  Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_x_th());
159 
160  TEUCHOS_ASSERT(r_!=Teuchos::null);
161 }
162 
163 // **********************************************************************
164 template<typename TRAITS,typename LO,typename GO>
165 void panzer::ScatterDirichletResidual_BlockedEpetra<panzer::Traits::Residual, TRAITS,LO,GO>::
166 evaluateFields(typename TRAITS::EvalData workset)
167 {
168  using Teuchos::RCP;
169  using Teuchos::ArrayRCP;
170  using Teuchos::ptrFromRef;
171  using Teuchos::rcp_dynamic_cast;
172 
173  using Thyra::VectorBase;
174  using Thyra::SpmdVectorBase;
175  using Thyra::ProductVectorBase;
176 
177  // for convenience pull out some objects from workset
178  std::string blockId = this->wda(workset).block_id;
179  const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
180 
181  // NOTE: A reordering of these loops will likely improve performance
182  // The "getGIDFieldOffsets may be expensive. However the
183  // "getElementGIDs" can be cheaper. However the lookup for LIDs
184  // may be more expensive!
185 
186  // loop over each field to be scattered
187  Teuchos::ArrayRCP<double> local_r;
188  Teuchos::ArrayRCP<double> local_dc;
189  for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
190  int rowIndexer = indexerIds_[fieldIndex];
191  int subFieldNum = subFieldIds_[fieldIndex];
192 
193  rcp_dynamic_cast<SpmdVectorBase<double> >(dirichletCounter_->getNonconstVectorBlock(rowIndexer))
194  ->getNonconstLocalData(ptrFromRef(local_dc));
195 
196  // grab local data for inputing
197  rcp_dynamic_cast<SpmdVectorBase<double> >(r_->getNonconstVectorBlock(rowIndexer))
198  ->getNonconstLocalData(ptrFromRef(local_r));
199 
200  auto subRowIndexer = rowIndexers_[rowIndexer];
201  auto LIDs = subRowIndexer->getLIDs();
202  auto LIDs_h = Kokkos::create_mirror_view(LIDs);
203  Kokkos::deep_copy(LIDs_h, LIDs);
204 
205  auto field = scatterFields_[fieldIndex].get_view();
206  auto field_h = Kokkos::create_mirror_view(field);
207  Kokkos::deep_copy(field_h, field);
208 
209  BCFieldType::array_type::HostMirror applyBC_h;
210  if(checkApplyBC_){
211  auto applyBC = applyBC_[fieldIndex].get_static_view();
212  applyBC_h = Kokkos::create_mirror_view(applyBC);
213  Kokkos::deep_copy(applyBC_h, applyBC);
214  }
215 
216  // scatter operation for each cell in workset
217  for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
218  std::size_t cellLocalId = localCellIds[worksetCellIndex];
219 
220  if (!scatterIC_) {
221  // this call "should" get the right ordering according to the Intrepid2 basis
222  const std::pair<std::vector<int>,std::vector<int> > & indicePair
223  = subRowIndexer->getGIDFieldOffsets_closure(blockId,subFieldNum, side_subcell_dim_, local_side_id_);
224  const std::vector<int> & elmtOffset = indicePair.first;
225  const std::vector<int> & basisIdMap = indicePair.second;
226 
227  // loop over basis functions
228  for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
229  int offset = elmtOffset[basis];
230  int lid = LIDs_h(cellLocalId, offset);
231  if(lid<0) // not on this processor!
232  continue;
233 
234  int basisId = basisIdMap[basis];
235 
236  if (checkApplyBC_ and !applyBC_h(worksetCellIndex,basisId))
237  continue;
238 
239  local_r[lid] = field_h(worksetCellIndex,basisId);
240 
241  // record that you set a dirichlet condition
242  local_dc[lid] = 1.0;
243  }
244  } else {
245  // this call "should" get the right ordering according to the Intrepid2 basis
246  const std::vector<int> & elmtOffset = subRowIndexer->getGIDFieldOffsets(blockId,subFieldNum);
247 
248  // loop over basis functions
249  for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
250  int offset = elmtOffset[basis];
251  int lid = LIDs_h(cellLocalId, offset);
252  if(lid<0) // not on this processor!
253  continue;
254 
255  local_r[lid] = field_h(worksetCellIndex,basis);
256 
257  // record that you set a dirichlet condition
258  local_dc[lid] = 1.0;
259  }
260  }
261  }
262  }
263 }
264 
265 // **********************************************************************
266 // Specialization: Tangent
267 // **********************************************************************
268 
269 
270 template<typename TRAITS,typename LO,typename GO>
271 panzer::ScatterDirichletResidual_BlockedEpetra<panzer::Traits::Tangent, TRAITS,LO,GO>::
272 ScatterDirichletResidual_BlockedEpetra(const std::vector<Teuchos::RCP<const GlobalIndexer> > & rIndexers,
273  const std::vector<Teuchos::RCP<const GlobalIndexer> > & cIndexers,
274  const Teuchos::ParameterList& p,
275  bool /* useDiscreteAdjoint */)
276  : rowIndexers_(rIndexers)
277  , colIndexers_(cIndexers)
278  , globalDataKey_("Residual Scatter Container")
279 {
280  std::string scatterName = p.get<std::string>("Scatter Name");
281  scatterHolder_ =
282  Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
283 
284  // get names to be evaluated
285  const std::vector<std::string>& names =
286  *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
287 
288  // grab map from evaluated names to field names
289  fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
290 
291  // determine if we are scattering an initial condition
292  scatterIC_ = p.isParameter("Scatter Initial Condition") ? p.get<bool>("Scatter Initial Condition") : false;
293 
294  Teuchos::RCP<PHX::DataLayout> dl = (!scatterIC_) ?
295  p.get< Teuchos::RCP<panzer::PureBasis> >("Basis")->functional :
296  p.get< Teuchos::RCP<const panzer::PureBasis> >("Basis")->functional;
297  if (!scatterIC_) {
298  side_subcell_dim_ = p.get<int>("Side Subcell Dimension");
299  local_side_id_ = p.get<int>("Local Side ID");
300  }
301 
302  // build the vector of fields that this is dependent on
303  scatterFields_.resize(names.size());
304  for (std::size_t eq = 0; eq < names.size(); ++eq) {
305  scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
306 
307  // tell the field manager that we depend on this field
308  this->addDependentField(scatterFields_[eq]);
309  }
310 
311  checkApplyBC_ = p.isParameter("Check Apply BC") ? p.get<bool>("Check Apply BC") : false;
312  if (checkApplyBC_) {
313  applyBC_.resize(names.size());
314  for (std::size_t eq = 0; eq < names.size(); ++eq) {
315  applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string("APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
316  this->addDependentField(applyBC_[eq]);
317  }
318  }
319 
320  // this is what this evaluator provides
321  this->addEvaluatedField(*scatterHolder_);
322 
323  if (p.isType<std::string>("Global Data Key"))
324  globalDataKey_ = p.get<std::string>("Global Data Key");
325 
326  this->setName(scatterName+" Scatter Tangent");
327 }
328 
329 // **********************************************************************
330 template<typename TRAITS,typename LO,typename GO>
331 void panzer::ScatterDirichletResidual_BlockedEpetra<panzer::Traits::Tangent, TRAITS,LO,GO>::
332 postRegistrationSetup(typename TRAITS::SetupData /* d */,
333  PHX::FieldManager<TRAITS>& /* fm */)
334 {
335  indexerIds_.resize(scatterFields_.size());
336  subFieldIds_.resize(scatterFields_.size());
337 
338  // load required field numbers for fast use
339  for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
340  // get field ID from DOF manager
341  std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
342 
343  indexerIds_[fd] = getFieldBlock(fieldName,rowIndexers_);
344  subFieldIds_[fd] = rowIndexers_[indexerIds_[fd]]->getFieldNum(fieldName);
345  }
346 
347  // get the number of nodes (Should be renamed basis)
348  num_nodes = scatterFields_[0].extent(1);
349 }
350 
351 // **********************************************************************
352 template<typename TRAITS,typename LO,typename GO>
353 void panzer::ScatterDirichletResidual_BlockedEpetra<panzer::Traits::Tangent, TRAITS,LO,GO>::
354 preEvaluate(typename TRAITS::PreEvalData d)
355 {
356  typedef BlockedEpetraLinearObjContainer BLOC;
357  typedef BlockedEpetraLinearObjContainer ELOC;
358 
359  using Teuchos::rcp_dynamic_cast;
360  using Thyra::ProductVectorBase;
361 
362  // extract dirichlet counter from container
363  Teuchos::RCP<BLOC> blockContainer
364  = Teuchos::rcp_dynamic_cast<BLOC>(d.gedc->getDataObject("Dirichlet Counter"),true);
365 
366  dirichletCounter_ = Teuchos::rcp_dynamic_cast<Thyra::ProductVectorBase<double> >(blockContainer->get_f(),true);
367  TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
368 
369  // extract linear object container
370  Teuchos::RCP<const BLOC> blockedContainer = Teuchos::rcp_dynamic_cast<const BLOC>(d.gedc->getDataObject(globalDataKey_));
371  Teuchos::RCP<const ELOC> epetraContainer = Teuchos::rcp_dynamic_cast<const ELOC>(d.gedc->getDataObject(globalDataKey_));
372 
373  // if its blocked do this
374  if(blockedContainer!=Teuchos::null)
375  r_ = (!scatterIC_) ?
376  rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer->get_f(),true) :
377  rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer->get_x(),true);
378  else if(epetraContainer!=Teuchos::null) // if its straight up epetra do this
379  r_ = (!scatterIC_) ?
380  Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_f_th()) :
381  Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_x_th());
382 
383  TEUCHOS_ASSERT(r_!=Teuchos::null);
384 }
385 
386 // **********************************************************************
387 template<typename TRAITS,typename LO,typename GO>
388 void panzer::ScatterDirichletResidual_BlockedEpetra<panzer::Traits::Tangent, TRAITS,LO,GO>::
389 evaluateFields(typename TRAITS::EvalData workset)
390 {
391  TEUCHOS_ASSERT(false);
392 
393  using Teuchos::RCP;
394  using Teuchos::ArrayRCP;
395  using Teuchos::ptrFromRef;
396  using Teuchos::rcp_dynamic_cast;
397 
398  using Thyra::VectorBase;
399  using Thyra::SpmdVectorBase;
400  using Thyra::ProductVectorBase;
401 
402  // for convenience pull out some objects from workset
403  std::string blockId = this->wda(workset).block_id;
404  const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
405 
406  // NOTE: A reordering of these loops will likely improve performance
407  // The "getGIDFieldOffsets may be expensive. However the
408  // "getElementGIDs" can be cheaper. However the lookup for LIDs
409  // may be more expensive!
410 
411  // loop over each field to be scattered
412  Teuchos::ArrayRCP<double> local_r;
413  Teuchos::ArrayRCP<double> local_dc;
414  for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
415  int rowIndexer = indexerIds_[fieldIndex];
416  int subFieldNum = subFieldIds_[fieldIndex];
417 
418  rcp_dynamic_cast<SpmdVectorBase<double> >(dirichletCounter_->getNonconstVectorBlock(rowIndexer))
419  ->getNonconstLocalData(ptrFromRef(local_dc));
420 
421  // grab local data for inputing
422  rcp_dynamic_cast<SpmdVectorBase<double> >(r_->getNonconstVectorBlock(rowIndexer))
423  ->getNonconstLocalData(ptrFromRef(local_r));
424 
425  auto subRowIndexer = rowIndexers_[rowIndexer];
426 
427  auto LIDs = subRowIndexer->getLIDs();
428  auto LIDs_h = Kokkos::create_mirror_view(LIDs);
429  Kokkos::deep_copy(LIDs_h, LIDs);
430 
431  auto field = scatterFields_[fieldIndex].get_view();
432  auto field_h = Kokkos::create_mirror_view(field);
433  Kokkos::deep_copy(field_h, field);
434 
435  BCFieldType::array_type::HostMirror applyBC_h;
436  if(checkApplyBC_){
437  auto applyBC = applyBC_[fieldIndex].get_static_view();
438  applyBC_h = Kokkos::create_mirror_view(applyBC);
439  Kokkos::deep_copy(applyBC_h, applyBC);
440  }
441 
442  // scatter operation for each cell in workset
443  for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
444  std::size_t cellLocalId = localCellIds[worksetCellIndex];
445 
446  if (!scatterIC_) {
447  // this call "should" get the right ordering according to the Intrepid2 basis
448  const std::pair<std::vector<int>,std::vector<int> > & indicePair
449  = subRowIndexer->getGIDFieldOffsets_closure(blockId,subFieldNum, side_subcell_dim_, local_side_id_);
450  const std::vector<int> & elmtOffset = indicePair.first;
451  const std::vector<int> & basisIdMap = indicePair.second;
452 
453  // loop over basis functions
454  for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
455  int offset = elmtOffset[basis];
456  int lid = LIDs_h(cellLocalId, offset);
457  if(lid<0) // not on this processor!
458  continue;
459 
460  int basisId = basisIdMap[basis];
461 
462  if (checkApplyBC_ and !applyBC_h(worksetCellIndex,basisId))
463  continue;
464 
465  local_r[lid] = field_h(worksetCellIndex,basisId).val();
466 
467  // record that you set a dirichlet condition
468  local_dc[lid] = 1.0;
469  }
470  } else {
471  // this call "should" get the right ordering according to the Intrepid2 basis
472  const std::vector<int> & elmtOffset = subRowIndexer->getGIDFieldOffsets(blockId,subFieldNum);
473 
474  // loop over basis functions
475  for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
476  int offset = elmtOffset[basis];
477  int lid = LIDs_h(cellLocalId, offset);
478  if(lid<0) // not on this processor!
479  continue;
480 
481  local_r[lid] = field_h(worksetCellIndex,basis).val();
482 
483  // record that you set a dirichlet condition
484  local_dc[lid] = 1.0;
485  }
486  }
487  }
488  }
489 }
490 
491 // **********************************************************************
492 // Specialization: Jacobian
493 // **********************************************************************
494 
495 template<typename TRAITS,typename LO,typename GO>
496 panzer::ScatterDirichletResidual_BlockedEpetra<panzer::Traits::Jacobian, TRAITS,LO,GO>::
497 ScatterDirichletResidual_BlockedEpetra(const std::vector<Teuchos::RCP<const GlobalIndexer> > & rIndexers,
498  const std::vector<Teuchos::RCP<const GlobalIndexer> > & cIndexers,
499  const Teuchos::ParameterList& p,
500  bool /* useDiscreteAdjoint */)
501  : rowIndexers_(rIndexers)
502  , colIndexers_(cIndexers)
503  , globalDataKey_("Residual Scatter Container")
504 {
505  std::string scatterName = p.get<std::string>("Scatter Name");
506  scatterHolder_ =
507  Teuchos::rcp(new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(new PHX::MDALayout<Dummy>(0))));
508 
509  // get names to be evaluated
510  const std::vector<std::string>& names =
511  *(p.get< Teuchos::RCP< std::vector<std::string> > >("Dependent Names"));
512 
513  // grab map from evaluated names to field names
514  fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >("Dependent Map");
515 
516  Teuchos::RCP<PHX::DataLayout> dl =
517  p.get< Teuchos::RCP<panzer::PureBasis> >("Basis")->functional;
518 
519  side_subcell_dim_ = p.get<int>("Side Subcell Dimension");
520  local_side_id_ = p.get<int>("Local Side ID");
521 
522  // build the vector of fields that this is dependent on
523  scatterFields_.resize(names.size());
524  for (std::size_t eq = 0; eq < names.size(); ++eq) {
525  scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
526 
527  // tell the field manager that we depend on this field
528  this->addDependentField(scatterFields_[eq]);
529  }
530 
531  checkApplyBC_ = p.get<bool>("Check Apply BC");
532  if (checkApplyBC_) {
533  applyBC_.resize(names.size());
534  for (std::size_t eq = 0; eq < names.size(); ++eq) {
535  applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string("APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
536  this->addDependentField(applyBC_[eq]);
537  }
538  }
539 
540  // this is what this evaluator provides
541  this->addEvaluatedField(*scatterHolder_);
542 
543  if (p.isType<std::string>("Global Data Key"))
544  globalDataKey_ = p.get<std::string>("Global Data Key");
545 
546  if(colIndexers_.size()==0)
547  colIndexers_ = rowIndexers_;
548 
549  this->setName(scatterName+" Scatter Residual (Jacobian)");
550 }
551 
552 // **********************************************************************
553 template<typename TRAITS,typename LO,typename GO>
554 void panzer::ScatterDirichletResidual_BlockedEpetra<panzer::Traits::Jacobian, TRAITS,LO,GO>::
555 postRegistrationSetup(typename TRAITS::SetupData /* d */,
556  PHX::FieldManager<TRAITS>& /* fm */)
557 {
558  indexerIds_.resize(scatterFields_.size());
559  subFieldIds_.resize(scatterFields_.size());
560 
561  // load required field numbers for fast use
562  for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
563  // get field ID from DOF manager
564  std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
565 
566  indexerIds_[fd] = getFieldBlock(fieldName,rowIndexers_);
567  subFieldIds_[fd] = rowIndexers_[indexerIds_[fd]]->getFieldNum(fieldName);
568  }
569 
570  // get the number of nodes (Should be renamed basis)
571  num_nodes = scatterFields_[0].extent(1);
572  num_eq = scatterFields_.size();
573 }
574 
575 // **********************************************************************
576 template<typename TRAITS,typename LO,typename GO>
577 void panzer::ScatterDirichletResidual_BlockedEpetra<panzer::Traits::Jacobian, TRAITS,LO,GO>::
578 preEvaluate(typename TRAITS::PreEvalData d)
579 {
580  typedef BlockedEpetraLinearObjContainer BLOC;
581 
582  using Teuchos::rcp_dynamic_cast;
583 
584  // extract dirichlet counter from container
585  Teuchos::RCP<const BLOC> blockContainer
586  = rcp_dynamic_cast<const BLOC>(d.gedc->getDataObject("Dirichlet Counter"),true);
587 
588  dirichletCounter_ = rcp_dynamic_cast<Thyra::ProductVectorBase<double> >(blockContainer->get_f(),true);
589  TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
590 
591  // extract linear object container
592  blockContainer = rcp_dynamic_cast<const BLOC>(d.gedc->getDataObject(globalDataKey_),true);
593  TEUCHOS_ASSERT(!Teuchos::is_null(blockContainer));
594 
595  r_ = rcp_dynamic_cast<Thyra::ProductVectorBase<double> >(blockContainer->get_f());
596  Jac_ = rcp_dynamic_cast<Thyra::BlockedLinearOpBase<double> >(blockContainer->get_A());
597 }
598 
599 // **********************************************************************
600 template<typename TRAITS,typename LO,typename GO>
601 void panzer::ScatterDirichletResidual_BlockedEpetra<panzer::Traits::Jacobian, TRAITS,LO,GO>::
602 evaluateFields(typename TRAITS::EvalData workset)
603 {
604  using Teuchos::RCP;
605  using Teuchos::ArrayRCP;
606  using Teuchos::ptrFromRef;
607  using Teuchos::rcp_dynamic_cast;
608 
609  using Thyra::SpmdVectorBase;
610 
611  // for convenience pull out some objects from workset
612  std::string blockId = this->wda(workset).block_id;
613  const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
614 
615  int numFieldBlocks = Teuchos::as<int>(colIndexers_.size());
616 
617  std::vector<int> blockOffsets;
618  computeBlockOffsets(blockId,colIndexers_,blockOffsets);
619 
620  std::unordered_map<std::pair<int,int>,Teuchos::RCP<Epetra_CrsMatrix>,panzer::pair_hash> jacEpetraBlocks;
621 
622  // NOTE: A reordering of these loops will likely improve performance
623  // The "getGIDFieldOffsets may be expensive. However the
624  // "getElementGIDs" can be cheaper. However the lookup for LIDs
625  // may be more expensive!
626 
627  for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
628  int rowIndexer = indexerIds_[fieldIndex];
629  int subFieldNum = subFieldIds_[fieldIndex];
630 
631  // loop over each field to be scattered
632  Teuchos::ArrayRCP<double> local_dc;
633  rcp_dynamic_cast<SpmdVectorBase<double> >(dirichletCounter_->getNonconstVectorBlock(rowIndexer))
634  ->getNonconstLocalData(ptrFromRef(local_dc));
635 
636  // grab local data for inputing
637  Teuchos::ArrayRCP<double> local_r;
638  if(r_!=Teuchos::null)
639  rcp_dynamic_cast<SpmdVectorBase<double> >(r_->getNonconstVectorBlock(rowIndexer))
640  ->getNonconstLocalData(ptrFromRef(local_r));
641 
642  auto subRowIndexer = rowIndexers_[rowIndexer];
643  auto subIndicePair = subRowIndexer->getGIDFieldOffsets_closure(blockId,subFieldNum, side_subcell_dim_, local_side_id_);
644  const std::vector<int> & subElmtOffset = subIndicePair.first;
645  const std::vector<int> & subBasisIdMap = subIndicePair.second;
646 
647  auto rLIDs = subRowIndexer->getLIDs();
648  auto rLIDs_h = Kokkos::create_mirror_view(rLIDs);
649  Kokkos::deep_copy(rLIDs_h, rLIDs);
650 
651  auto field = scatterFields_[fieldIndex].get_view();
652  auto field_h = Kokkos::create_mirror_view(field);
653  Kokkos::deep_copy(field_h, field);
654 
655  BCFieldType::array_type::HostMirror applyBC_h;
656  if(checkApplyBC_){
657  auto applyBC = applyBC_[fieldIndex].get_static_view();
658  applyBC_h = Kokkos::create_mirror_view(applyBC);
659  Kokkos::deep_copy(applyBC_h, applyBC);
660  }
661 
662  // scatter operation for each cell in workset
663  for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
664  std::size_t cellLocalId = localCellIds[worksetCellIndex];
665 
666  // loop over basis functions
667  for(std::size_t basis=0;basis<subElmtOffset.size();basis++) {
668  int offset = subElmtOffset[basis];
669  int lid = rLIDs_h(cellLocalId, offset);
670  if(lid<0) // not on this processor
671  continue;
672 
673  int basisId = subBasisIdMap[basis];
674 
675  if (checkApplyBC_ and !applyBC_h(worksetCellIndex,basisId))
676  continue;
677 
678  // zero out matrix row
679  for(int colIndexer=0;colIndexer<numFieldBlocks;colIndexer++) {
680  int start = blockOffsets[colIndexer];
681  int end = blockOffsets[colIndexer+1];
682 
683  if(end-start<=0)
684  continue;
685 
686  // check hash table for jacobian sub block
687  std::pair<int,int> blockIndex = std::make_pair(rowIndexer,colIndexer);
688  Teuchos::RCP<Epetra_CrsMatrix> subJac = jacEpetraBlocks[blockIndex];
689  // if you didn't find one before, add it to the hash table
690  if(subJac==Teuchos::null) {
691  Teuchos::RCP<Thyra::LinearOpBase<double> > tOp = Jac_->getNonconstBlock(blockIndex.first,blockIndex.second);
692 
693  // block operator is null, don't do anything (it is excluded)
694  if(Teuchos::is_null(tOp))
695  continue;
696 
697  Teuchos::RCP<Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*tOp);
698  subJac = rcp_dynamic_cast<Epetra_CrsMatrix>(eOp,true);
699  jacEpetraBlocks[blockIndex] = subJac;
700  }
701 
702  int numEntries = 0;
703  int * rowIndices = 0;
704  double * rowValues = 0;
705 
706  subJac->ExtractMyRowView(lid,numEntries,rowValues,rowIndices);
707 
708  for(int i=0;i<numEntries;i++)
709  rowValues[i] = 0.0;
710  }
711 
712  const ScalarT scatterField = field_h(worksetCellIndex,basisId);
713 
714  if(r_!=Teuchos::null)
715  local_r[lid] = scatterField.val();
716  local_dc[lid] = 1.0; // mark row as dirichlet
717 
718  // loop over the sensitivity indices: all DOFs on a cell
719  std::vector<double> jacRow(scatterField.size(),0.0);
720 
721  for(int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
722  jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex);
723 
724  for(int colIndexer=0;colIndexer<numFieldBlocks;colIndexer++) {
725  int start = blockOffsets[colIndexer];
726  int end = blockOffsets[colIndexer+1];
727 
728  if(end-start<=0)
729  continue;
730 
731  auto subColIndexer = colIndexers_[colIndexer];
732 
733  auto cLIDs = subColIndexer->getElementLIDs(cellLocalId);
734  auto cLIDs_h = Kokkos::create_mirror_view(cLIDs);
735  Kokkos::deep_copy(cLIDs_h, cLIDs);
736 
737  TEUCHOS_ASSERT(end-start==Teuchos::as<int>(cLIDs.size()));
738 
739  // check hash table for jacobian sub block
740  std::pair<int,int> blockIndex = std::make_pair(rowIndexer,colIndexer);
741  Teuchos::RCP<Epetra_CrsMatrix> subJac = jacEpetraBlocks[blockIndex];
742 
743  // if you didn't find one before, add it to the hash table
744  if(subJac==Teuchos::null) {
745  Teuchos::RCP<Thyra::LinearOpBase<double> > tOp = Jac_->getNonconstBlock(blockIndex.first,blockIndex.second);
746 
747  // block operator is null, don't do anything (it is excluded)
748  if(Teuchos::is_null(tOp))
749  continue;
750 
751  Teuchos::RCP<Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*tOp);
752  subJac = rcp_dynamic_cast<Epetra_CrsMatrix>(eOp,true);
753  jacEpetraBlocks[blockIndex] = subJac;
754  }
755 
756  // Sum Jacobian
757  int err = subJac->ReplaceMyValues(lid, end-start, &jacRow[start],&cLIDs_h[0]);
758  if(err!=0) {
759  std::stringstream ss;
760  ss << "Failed inserting row: " << " (" << lid << "): ";
761  for(int i=0;i<end-start;i++)
762  ss << cLIDs_h[i] << " ";
763  ss << std::endl;
764  ss << "Into block " << rowIndexer << ", " << colIndexer << std::endl;
765 
766  ss << "scatter field = ";
767  scatterFields_[fieldIndex].print(ss);
768  ss << std::endl;
769 
770  TEUCHOS_TEST_FOR_EXCEPTION(err!=0,std::runtime_error,ss.str());
771  }
772 
773  }
774  }
775  }
776  }
777 }
778 
779 // **********************************************************************
780 
781 #endif
panzer::pair_hash
Definition: Panzer_HashUtils.hpp:46
Teuchos::ParameterList::get
T & get(const std::string &name, T def_value)
Teuchos::is_null
bool is_null(const std::shared_ptr< T > &p)
Panzer_PureBasis.hpp
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
panzer::ScatterDirichletResidual_BlockedEpetra::postRegistrationSetup
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &vm)
Definition: Panzer_ScatterDirichletResidual_BlockedEpetra_decl.hpp:57
panzer::computeBlockOffsets
void computeBlockOffsets(const std::string &blockId, const std::vector< Teuchos::RCP< GlobalIndexer >> &ugis, std::vector< int > &blockOffsets)
Definition: Panzer_GlobalIndexer_Utilities.cpp:54
Epetra_CrsMatrix::ExtractMyRowView
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
Teuchos::ParameterList::isParameter
bool isParameter(const std::string &name) const
Panzer_BlockedEpetraLinearObjContainer.hpp
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
panzer::ScatterDirichletResidual_BlockedEpetra::evaluateFields
void evaluateFields(typename TRAITS::EvalData d)
Definition: Panzer_ScatterDirichletResidual_BlockedEpetra_decl.hpp:59
panzer::getFieldBlock
int getFieldBlock(const std::string &fieldName, const std::vector< Teuchos::RCP< const GlobalIndexer >> &ugis)
Definition: Panzer_GlobalIndexer_Utilities.cpp:28
panzer::ScatterDirichletResidual_BlockedEpetra::ScatterDirichletResidual_BlockedEpetra
ScatterDirichletResidual_BlockedEpetra(const Teuchos::ParameterList &p)
Definition: Panzer_ScatterDirichletResidual_BlockedEpetra_decl.hpp:52
Panzer_GlobalIndexer_Utilities.hpp
PHX::Evaluator::preEvaluate
virtual void preEvaluate(typename Traits::PreEvalData d)=0
panzer::ScatterDirichletResidual_BlockedEpetra< panzer::Traits::Jacobian, TRAITS, LO, GO >::ScalarT
panzer::Traits::Jacobian::ScalarT ScalarT
Definition: Panzer_ScatterDirichletResidual_BlockedEpetra_decl.hpp:249
field
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...
Definition: Panzer_Integrator_CurlBasisDotVector_impl.hpp:559
Epetra_CrsMatrix
Panzer_BlockedDOFManager.hpp
Teuchos::ParameterList::isType
bool isType(const std::string &name) const
Panzer_GlobalEvaluationDataContainer.hpp
Panzer_GlobalIndexer.hpp
TEUCHOS_ASSERT
#define TEUCHOS_ASSERT(assertion_test)
Epetra_CrsMatrix::ReplaceMyValues
int ReplaceMyValues(int MyRow, int NumEntries, const double *Values, const int *Indices)
end
TransListIter end
Teuchos_RCP.hpp
Generated on Fri Dec 20 2024 09:16:52 for Panzer by   doxygen 1.8.5