43 #ifndef PANZER_SCATTER_DIRICHLET_RESIDUAL_EPETRA_IMPL_HPP 44 #define PANZER_SCATTER_DIRICHLET_RESIDUAL_EPETRA_IMPL_HPP 46 #include "Teuchos_RCP.hpp" 47 #include "Teuchos_Assert.hpp" 49 #include "Phalanx_DataLayout.hpp" 51 #include "Epetra_Map.h" 52 #include "Epetra_Vector.h" 53 #include "Epetra_CrsMatrix.h" 62 #include "Phalanx_DataLayout_MDALayout.hpp" 64 #include "Teuchos_FancyOStream.hpp" 71 template<
typename TRAITS,
typename LO,
typename GO>
74 const Teuchos::RCP<const panzer::GlobalIndexer> & ,
75 const Teuchos::ParameterList& p)
76 : globalIndexer_(indexer)
77 , globalDataKey_(
"Residual Scatter Container")
79 std::string scatterName = p.get<std::string>(
"Scatter Name");
81 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
84 const std::vector<std::string>& names =
85 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
88 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
91 scatterIC_ = p.isParameter(
"Scatter Initial Condition") ? p.get<
bool>(
"Scatter Initial Condition") :
false;
93 Teuchos::RCP<PHX::DataLayout> dl = (!scatterIC_) ?
94 p.get< Teuchos::RCP<panzer::PureBasis> >(
"Basis")->functional :
95 p.get< Teuchos::RCP<const panzer::PureBasis> >(
"Basis")->functional;
97 side_subcell_dim_ = p.get<
int>(
"Side Subcell Dimension");
98 local_side_id_ = p.get<
int>(
"Local Side ID");
102 scatterFields_.resize(names.size());
103 for (std::size_t eq = 0; eq < names.size(); ++eq) {
104 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
107 this->addDependentField(scatterFields_[eq]);
110 checkApplyBC_ = p.isParameter(
"Check Apply BC") ? p.get<
bool>(
"Check Apply BC") :
false;
112 applyBC_.resize(names.size());
113 for (std::size_t eq = 0; eq < names.size(); ++eq) {
114 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string(
"APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
115 this->addDependentField(applyBC_[eq]);
120 this->addEvaluatedField(*scatterHolder_);
122 if (p.isType<std::string>(
"Global Data Key"))
123 globalDataKey_ = p.get<std::string>(
"Global Data Key");
125 this->setName(scatterName+
" Scatter Residual");
129 template<
typename TRAITS,
typename LO,
typename GO>
134 fieldIds_.resize(scatterFields_.size());
137 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
139 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
140 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
144 num_nodes = scatterFields_[0].extent(1);
148 template<
typename TRAITS,
typename LO,
typename GO>
155 if(epetraContainer_==Teuchos::null) {
157 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<
LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),
true)->getGhostedLOC();
160 dirichletCounter_ = Teuchos::null;
164 Teuchos::RCP<EpetraLinearObjContainer> epetraContainer
167 dirichletCounter_ = epetraContainer->
get_f();
168 TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
173 template<
typename TRAITS,
typename LO,
typename GO>
178 std::string blockId = this->wda(workset).block_id;
179 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
181 Teuchos::RCP<const EpetraLinearObjContainer> epetraContainer = epetraContainer_;
182 Teuchos::RCP<Epetra_Vector> r = (!scatterIC_) ?
183 epetraContainer->get_f() :
184 epetraContainer->get_x();
191 auto LIDs = globalIndexer_->getLIDs();
192 auto LIDs_h = Kokkos::create_mirror_view(LIDs);
193 Kokkos::deep_copy(LIDs_h, LIDs);
197 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
198 int fieldNum = fieldIds_[fieldIndex];
199 auto field = PHX::as_view(scatterFields_[fieldIndex]);
200 auto field_h = Kokkos::create_mirror_view(
field);
201 Kokkos::deep_copy(field_h,
field);
203 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
204 std::size_t cellLocalId = localCellIds[worksetCellIndex];
208 const std::pair<std::vector<int>,std::vector<int> > & indicePair
209 = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldNum, side_subcell_dim_, local_side_id_);
210 const std::vector<int> & elmtOffset = indicePair.first;
211 const std::vector<int> & basisIdMap = indicePair.second;
214 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
215 int offset = elmtOffset[basis];
216 int lid = LIDs_h(cellLocalId, offset);
220 int basisId = basisIdMap[basis];
223 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
226 (*r)[lid] = field_h(worksetCellIndex,basisId);
229 if(dirichletCounter_!=Teuchos::null)
230 (*dirichletCounter_)[lid] = 1.0;
234 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
237 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
238 int offset = elmtOffset[basis];
239 int lid = LIDs_h(cellLocalId, offset);
243 (*r)[lid] = field_h(worksetCellIndex,basis);
246 if(dirichletCounter_!=Teuchos::null)
247 (*dirichletCounter_)[lid] = 1.0;
259 template<
typename TRAITS,
typename LO,
typename GO>
262 const Teuchos::RCP<const panzer::GlobalIndexer> & ,
263 const Teuchos::ParameterList& p)
264 : globalIndexer_(indexer)
265 , globalDataKey_(
"Residual Scatter Container")
267 std::string scatterName = p.get<std::string>(
"Scatter Name");
269 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
272 const std::vector<std::string>& names =
273 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
276 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
279 scatterIC_ = p.isParameter(
"Scatter Initial Condition") ? p.get<
bool>(
"Scatter Initial Condition") :
false;
281 Teuchos::RCP<PHX::DataLayout> dl = (!scatterIC_) ?
282 p.get< Teuchos::RCP<panzer::PureBasis> >(
"Basis")->functional :
283 p.get< Teuchos::RCP<const panzer::PureBasis> >(
"Basis")->functional;
285 side_subcell_dim_ = p.get<
int>(
"Side Subcell Dimension");
286 local_side_id_ = p.get<
int>(
"Local Side ID");
290 scatterFields_.resize(names.size());
291 for (std::size_t eq = 0; eq < names.size(); ++eq) {
292 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
295 this->addDependentField(scatterFields_[eq]);
298 checkApplyBC_ = p.isParameter(
"Check Apply BC") ? p.get<
bool>(
"Check Apply BC") :
false;
300 applyBC_.resize(names.size());
301 for (std::size_t eq = 0; eq < names.size(); ++eq) {
302 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string(
"APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
303 this->addDependentField(applyBC_[eq]);
308 this->addEvaluatedField(*scatterHolder_);
310 if (p.isType<std::string>(
"Global Data Key"))
311 globalDataKey_ = p.get<std::string>(
"Global Data Key");
313 this->setName(scatterName+
" Scatter Tangent");
317 template<
typename TRAITS,
typename LO,
typename GO>
322 fieldIds_.resize(scatterFields_.size());
325 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
327 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
328 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
332 num_nodes = scatterFields_[0].extent(1);
336 template<
typename TRAITS,
typename LO,
typename GO>
343 if(epetraContainer_==Teuchos::null) {
345 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<
LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),
true)->getGhostedLOC();
348 dirichletCounter_ = Teuchos::null;
352 Teuchos::RCP<EpetraLinearObjContainer> epetraContainer
355 dirichletCounter_ = epetraContainer->
get_f();
356 TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
360 using Teuchos::rcp_dynamic_cast;
363 std::vector<std::string> activeParameters =
367 TEUCHOS_ASSERT(!scatterIC_);
368 dfdp_vectors_.clear();
369 for(std::size_t i=0;i<activeParameters.size();i++) {
370 RCP<Epetra_Vector> vec = rcp_dynamic_cast<
EpetraLinearObjContainer>(d.gedc->getDataObject(activeParameters[i]),
true)->get_f();
371 dfdp_vectors_.push_back(vec);
376 template<
typename TRAITS,
typename LO,
typename GO>
381 std::string blockId = this->wda(workset).block_id;
382 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
384 Teuchos::RCP<const EpetraLinearObjContainer> epetraContainer = epetraContainer_;
385 Teuchos::RCP<Epetra_Vector> r = (!scatterIC_) ?
386 epetraContainer->get_f() :
387 epetraContainer->get_x();
394 auto LIDs = globalIndexer_->getLIDs();
395 auto LIDs_h = Kokkos::create_mirror_view(LIDs);
396 Kokkos::deep_copy(LIDs_h, LIDs);
397 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
398 int fieldNum = fieldIds_[fieldIndex];
399 auto scatterField_h = Kokkos::create_mirror_view(scatterFields_[fieldIndex].get_static_view());
400 Kokkos::deep_copy(scatterField_h, scatterFields_[fieldIndex].get_static_view());
403 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
404 std::size_t cellLocalId = localCellIds[worksetCellIndex];
408 const std::pair<std::vector<int>,std::vector<int> > & indicePair
409 = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldNum, side_subcell_dim_, local_side_id_);
410 const std::vector<int> & elmtOffset = indicePair.first;
411 const std::vector<int> & basisIdMap = indicePair.second;
414 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
415 int offset = elmtOffset[basis];
416 int lid = LIDs_h(cellLocalId, offset);
420 int basisId = basisIdMap[basis];
423 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
426 ScalarT value = scatterField_h(worksetCellIndex,basisId);
431 for(std::size_t d=0;d<dfdp_vectors_.size();d++)
432 (*dfdp_vectors_[d])[lid] = 0.0;
434 for(
int d=0;d<value.size();d++) {
435 (*dfdp_vectors_[d])[lid] = value.fastAccessDx(d);
439 if(dirichletCounter_!=Teuchos::null)
440 (*dirichletCounter_)[lid] = 1.0;
444 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
447 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
448 int offset = elmtOffset[basis];
449 int lid = LIDs_h(cellLocalId, offset);
453 ScalarT value = scatterField_h(worksetCellIndex,basis);
458 for(std::size_t d=0;d<dfdp_vectors_.size();d++)
459 (*dfdp_vectors_[d])[lid] = 0.0;
461 for(
int d=0;d<value.size();d++) {
462 (*dfdp_vectors_[d])[lid] = value.fastAccessDx(d);
466 if(dirichletCounter_!=Teuchos::null)
467 (*dirichletCounter_)[lid] = 1.0;
478 template<
typename TRAITS,
typename LO,
typename GO>
481 const Teuchos::RCP<const panzer::GlobalIndexer> & cIndexer,
482 const Teuchos::ParameterList& p)
483 : globalIndexer_(indexer)
484 , colGlobalIndexer_(cIndexer)
485 , globalDataKey_(
"Residual Scatter Container")
486 , preserveDiagonal_(false)
488 std::string scatterName = p.get<std::string>(
"Scatter Name");
490 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
493 const std::vector<std::string>& names =
494 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
497 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
499 Teuchos::RCP<PHX::DataLayout> dl =
500 p.get< Teuchos::RCP<panzer::PureBasis> >(
"Basis")->functional;
502 side_subcell_dim_ = p.get<
int>(
"Side Subcell Dimension");
503 local_side_id_ = p.get<
int>(
"Local Side ID");
506 scatterFields_.resize(names.size());
507 for (std::size_t eq = 0; eq < names.size(); ++eq) {
508 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
511 this->addDependentField(scatterFields_[eq]);
514 checkApplyBC_ = p.get<
bool>(
"Check Apply BC");
516 applyBC_.resize(names.size());
517 for (std::size_t eq = 0; eq < names.size(); ++eq) {
518 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string(
"APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
519 this->addDependentField(applyBC_[eq]);
524 this->addEvaluatedField(*scatterHolder_);
526 if (p.isType<std::string>(
"Global Data Key"))
527 globalDataKey_ = p.get<std::string>(
"Global Data Key");
529 if (p.isType<
bool>(
"Preserve Diagonal"))
530 preserveDiagonal_ = p.get<
bool>(
"Preserve Diagonal");
532 this->setName(scatterName+
" Scatter Residual (Jacobian)");
536 template<
typename TRAITS,
typename LO,
typename GO>
541 fieldIds_.resize(scatterFields_.size());
544 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
546 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
547 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
551 num_nodes = scatterFields_[0].extent(1);
552 num_eq = scatterFields_.size();
556 template<
typename TRAITS,
typename LO,
typename GO>
563 if(epetraContainer_==Teuchos::null) {
565 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<
LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),
true)->getGhostedLOC();
568 dirichletCounter_ = Teuchos::null;
572 Teuchos::RCP<GlobalEvaluationData> dataContainer = d.gedc->getDataObject(
"Dirichlet Counter");
573 Teuchos::RCP<EpetraLinearObjContainer> epetraContainer = Teuchos::rcp_dynamic_cast<
EpetraLinearObjContainer>(dataContainer,
true);
575 dirichletCounter_ = epetraContainer->
get_f();
576 TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
581 template<
typename TRAITS,
typename LO,
typename GO>
585 Kokkos::View<const int*, Kokkos::LayoutRight, PHX::Device> cLIDs, rLIDs;
587 bool useColumnIndexer = colGlobalIndexer_!=Teuchos::null;
588 const Teuchos::RCP<const panzer::GlobalIndexer>&
589 colGlobalIndexer = useColumnIndexer ? colGlobalIndexer_ : globalIndexer_;
592 std::string blockId = this->wda(workset).block_id;
593 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
595 Teuchos::RCP<const EpetraLinearObjContainer> epetraContainer = epetraContainer_;
596 TEUCHOS_ASSERT(epetraContainer!=Teuchos::null);
597 Teuchos::RCP<Epetra_Vector> r = epetraContainer->get_f();
598 Teuchos::RCP<Epetra_CrsMatrix> Jac = epetraContainer->get_A();
606 auto LIDs = globalIndexer_->getLIDs();
607 auto colLIDs = colGlobalIndexer->getLIDs();
608 auto LIDs_h = Kokkos::create_mirror_view(LIDs);
609 auto colLIDs_h = Kokkos::create_mirror_view(colLIDs);
610 Kokkos::deep_copy(LIDs_h, LIDs);
611 Kokkos::deep_copy(colLIDs_h, colLIDs);
613 std::vector<
typename decltype(scatterFields_[0].get_static_view())::HostMirror> scatterFields_h;
614 for ( std::size_t i=0; i< scatterFields_.size(); ++i) {
615 scatterFields_h.push_back(Kokkos::create_mirror_view(scatterFields_[i].get_static_view()));
616 Kokkos::deep_copy(scatterFields_h[i], scatterFields_[i].get_static_view());
619 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
620 std::size_t cellLocalId = localCellIds[worksetCellIndex];
622 gidCount = colGlobalIndexer->getElementBlockGIDCount(blockId);
625 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
626 int fieldNum = fieldIds_[fieldIndex];
629 const std::pair<std::vector<int>,std::vector<int> > & indicePair
630 = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldNum, side_subcell_dim_, local_side_id_);
631 const std::vector<int> & elmtOffset = indicePair.first;
632 const std::vector<int> & basisIdMap = indicePair.second;
635 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
636 int offset = elmtOffset[basis];
637 int row = LIDs_h(cellLocalId, offset);
641 int basisId = basisIdMap[basis];
644 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
650 int * rowIndices = 0;
651 double * rowValues = 0;
653 Jac->ExtractMyRowView(row,numEntries,rowValues,rowIndices);
655 for(
int i=0;i<numEntries;i++) {
656 if(preserveDiagonal_) {
657 if(row!=rowIndices[i])
666 const ScalarT scatterField = (scatterFields_h[fieldIndex])(worksetCellIndex,basisId);
669 (*r)[row] = scatterField.val();
670 if(dirichletCounter_!=Teuchos::null) {
672 (*dirichletCounter_)[row] = 1.0;
676 std::vector<double> jacRow(scatterField.size(),0.0);
678 if(!preserveDiagonal_) {
679 int err = Jac->ReplaceMyValues(row, gidCount, scatterField.dx(),
680 &colLIDs_h(cellLocalId,0));
681 TEUCHOS_ASSERT(err==0);
panzer::Traits::Tangent::ScalarT ScalarT
panzer::Traits::Jacobian::ScalarT ScalarT
PHX::MDField< ScalarT, panzer::Cell, panzer::BASIS > field
A field to which we'll contribute, or in which we'll store, the result of computing this integral...
const Teuchos::RCP< Epetra_Vector > get_f() const
Pushes residual values into the residual vector for a Newton-based solve.