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>
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 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
192 std::size_t cellLocalId = localCellIds[worksetCellIndex];
194 auto LIDs = globalIndexer_->getElementLIDs(cellLocalId);
197 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
198 int fieldNum = fieldIds_[fieldIndex];
202 const std::pair<std::vector<int>,std::vector<int> > & indicePair
203 = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldNum, side_subcell_dim_, local_side_id_);
204 const std::vector<int> & elmtOffset = indicePair.first;
205 const std::vector<int> & basisIdMap = indicePair.second;
208 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
209 int offset = elmtOffset[basis];
210 int lid = LIDs[offset];
214 int basisId = basisIdMap[basis];
217 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
220 (*r)[lid] = (scatterFields_[fieldIndex])(worksetCellIndex,basisId);
223 if(dirichletCounter_!=Teuchos::null)
224 (*dirichletCounter_)[lid] = 1.0;
228 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
231 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
232 int offset = elmtOffset[basis];
233 int lid = LIDs[offset];
237 (*r)[lid] = (scatterFields_[fieldIndex])(worksetCellIndex,basis);
240 if(dirichletCounter_!=Teuchos::null)
241 (*dirichletCounter_)[lid] = 1.0;
253 template<
typename TRAITS,
typename LO,
typename GO>
257 const Teuchos::ParameterList& p)
258 : globalIndexer_(indexer)
259 , globalDataKey_(
"Residual Scatter Container")
261 std::string scatterName = p.get<std::string>(
"Scatter Name");
263 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
266 const std::vector<std::string>& names =
267 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
270 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
273 scatterIC_ = p.isParameter(
"Scatter Initial Condition") ? p.get<
bool>(
"Scatter Initial Condition") :
false;
275 Teuchos::RCP<PHX::DataLayout> dl = (!scatterIC_) ?
276 p.get< Teuchos::RCP<panzer::PureBasis> >(
"Basis")->functional :
277 p.get< Teuchos::RCP<const panzer::PureBasis> >(
"Basis")->functional;
279 side_subcell_dim_ = p.get<
int>(
"Side Subcell Dimension");
280 local_side_id_ = p.get<
int>(
"Local Side ID");
284 scatterFields_.resize(names.size());
285 for (std::size_t eq = 0; eq < names.size(); ++eq) {
286 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
289 this->addDependentField(scatterFields_[eq]);
292 checkApplyBC_ = p.isParameter(
"Check Apply BC") ? p.get<
bool>(
"Check Apply BC") :
false;
294 applyBC_.resize(names.size());
295 for (std::size_t eq = 0; eq < names.size(); ++eq) {
296 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string(
"APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
297 this->addDependentField(applyBC_[eq]);
302 this->addEvaluatedField(*scatterHolder_);
304 if (p.isType<std::string>(
"Global Data Key"))
305 globalDataKey_ = p.get<std::string>(
"Global Data Key");
307 this->setName(scatterName+
" Scatter Tangent");
311 template<
typename TRAITS,
typename LO,
typename GO>
316 fieldIds_.resize(scatterFields_.size());
319 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
321 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
322 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
326 num_nodes = scatterFields_[0].extent(1);
330 template<
typename TRAITS,
typename LO,
typename GO>
337 if(epetraContainer_==Teuchos::null) {
339 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<
LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),
true)->getGhostedLOC();
342 dirichletCounter_ = Teuchos::null;
346 Teuchos::RCP<EpetraLinearObjContainer> epetraContainer
349 dirichletCounter_ = epetraContainer->
get_f();
350 TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
354 using Teuchos::rcp_dynamic_cast;
357 std::vector<std::string> activeParameters =
361 TEUCHOS_ASSERT(!scatterIC_);
362 dfdp_vectors_.clear();
363 for(std::size_t i=0;i<activeParameters.size();i++) {
364 RCP<Epetra_Vector> vec = rcp_dynamic_cast<
EpetraLinearObjContainer>(d.gedc->getDataObject(activeParameters[i]),
true)->get_f();
365 dfdp_vectors_.push_back(vec);
370 template<
typename TRAITS,
typename LO,
typename GO>
375 std::string blockId = this->wda(workset).block_id;
376 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
378 Teuchos::RCP<const EpetraLinearObjContainer> epetraContainer = epetraContainer_;
379 Teuchos::RCP<Epetra_Vector> r = (!scatterIC_) ?
380 epetraContainer->get_f() :
381 epetraContainer->get_x();
388 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
389 std::size_t cellLocalId = localCellIds[worksetCellIndex];
392 auto LIDs = globalIndexer_->getElementLIDs(cellLocalId);
395 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
396 int fieldNum = fieldIds_[fieldIndex];
400 const std::pair<std::vector<int>,std::vector<int> > & indicePair
401 = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldNum, side_subcell_dim_, local_side_id_);
402 const std::vector<int> & elmtOffset = indicePair.first;
403 const std::vector<int> & basisIdMap = indicePair.second;
406 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
407 int offset = elmtOffset[basis];
408 int lid = LIDs[offset];
412 int basisId = basisIdMap[basis];
415 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
418 ScalarT value = (scatterFields_[fieldIndex])(worksetCellIndex,basisId);
423 for(std::size_t d=0;d<dfdp_vectors_.size();d++)
424 (*dfdp_vectors_[d])[lid] = 0.0;
426 for(
int d=0;d<value.size();d++) {
427 (*dfdp_vectors_[d])[lid] = value.fastAccessDx(d);
431 if(dirichletCounter_!=Teuchos::null)
432 (*dirichletCounter_)[lid] = 1.0;
436 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
439 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
440 int offset = elmtOffset[basis];
441 int lid = LIDs[offset];
445 ScalarT value = (scatterFields_[fieldIndex])(worksetCellIndex,basis);
450 for(std::size_t d=0;d<dfdp_vectors_.size();d++)
451 (*dfdp_vectors_[d])[lid] = 0.0;
453 for(
int d=0;d<value.size();d++) {
454 (*dfdp_vectors_[d])[lid] = value.fastAccessDx(d);
458 if(dirichletCounter_!=Teuchos::null)
459 (*dirichletCounter_)[lid] = 1.0;
470 template<
typename TRAITS,
typename LO,
typename GO>
474 const Teuchos::ParameterList& p)
475 : globalIndexer_(indexer)
476 , colGlobalIndexer_(cIndexer)
477 , globalDataKey_(
"Residual Scatter Container")
478 , preserveDiagonal_(false)
480 std::string scatterName = p.get<std::string>(
"Scatter Name");
482 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
485 const std::vector<std::string>& names =
486 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
489 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
491 Teuchos::RCP<PHX::DataLayout> dl =
492 p.get< Teuchos::RCP<panzer::PureBasis> >(
"Basis")->functional;
494 side_subcell_dim_ = p.get<
int>(
"Side Subcell Dimension");
495 local_side_id_ = p.get<
int>(
"Local Side ID");
498 scatterFields_.resize(names.size());
499 for (std::size_t eq = 0; eq < names.size(); ++eq) {
500 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
503 this->addDependentField(scatterFields_[eq]);
506 checkApplyBC_ = p.get<
bool>(
"Check Apply BC");
508 applyBC_.resize(names.size());
509 for (std::size_t eq = 0; eq < names.size(); ++eq) {
510 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string(
"APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
511 this->addDependentField(applyBC_[eq]);
516 this->addEvaluatedField(*scatterHolder_);
518 if (p.isType<std::string>(
"Global Data Key"))
519 globalDataKey_ = p.get<std::string>(
"Global Data Key");
521 if (p.isType<
bool>(
"Preserve Diagonal"))
522 preserveDiagonal_ = p.get<
bool>(
"Preserve Diagonal");
524 this->setName(scatterName+
" Scatter Residual (Jacobian)");
528 template<
typename TRAITS,
typename LO,
typename GO>
533 fieldIds_.resize(scatterFields_.size());
536 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
538 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
539 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
543 num_nodes = scatterFields_[0].extent(1);
544 num_eq = scatterFields_.size();
548 template<
typename TRAITS,
typename LO,
typename GO>
555 if(epetraContainer_==Teuchos::null) {
557 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<
LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),
true)->getGhostedLOC();
560 dirichletCounter_ = Teuchos::null;
564 Teuchos::RCP<GlobalEvaluationData> dataContainer = d.gedc->getDataObject(
"Dirichlet Counter");
565 Teuchos::RCP<EpetraLinearObjContainer> epetraContainer = Teuchos::rcp_dynamic_cast<
EpetraLinearObjContainer>(dataContainer,
true);
567 dirichletCounter_ = epetraContainer->
get_f();
568 TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
573 template<
typename TRAITS,
typename LO,
typename GO>
577 Kokkos::View<const int*, Kokkos::LayoutRight, PHX::Device> cLIDs, rLIDs;
579 bool useColumnIndexer = colGlobalIndexer_!=Teuchos::null;
582 std::string blockId = this->wda(workset).block_id;
583 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
585 Teuchos::RCP<const EpetraLinearObjContainer> epetraContainer = epetraContainer_;
586 TEUCHOS_ASSERT(epetraContainer!=Teuchos::null);
587 Teuchos::RCP<Epetra_Vector> r = epetraContainer->get_f();
588 Teuchos::RCP<Epetra_CrsMatrix> Jac = epetraContainer->get_A();
595 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
596 std::size_t cellLocalId = localCellIds[worksetCellIndex];
598 rLIDs = globalIndexer_->getElementLIDs(cellLocalId);
599 gidCount = globalIndexer_->getElementBlockGIDCount(blockId);
602 cLIDs = colGlobalIndexer_->getElementLIDs(cellLocalId);
603 gidCount = colGlobalIndexer_->getElementBlockGIDCount(blockId);
609 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
610 int fieldNum = fieldIds_[fieldIndex];
613 const std::pair<std::vector<int>,std::vector<int> > & indicePair
614 = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldNum, side_subcell_dim_, local_side_id_);
615 const std::vector<int> & elmtOffset = indicePair.first;
616 const std::vector<int> & basisIdMap = indicePair.second;
619 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
620 int offset = elmtOffset[basis];
621 int row = rLIDs[offset];
625 int basisId = basisIdMap[basis];
628 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
634 int * rowIndices = 0;
635 double * rowValues = 0;
637 Jac->ExtractMyRowView(row,numEntries,rowValues,rowIndices);
639 for(
int i=0;i<numEntries;i++) {
640 if(preserveDiagonal_) {
641 if(row!=rowIndices[i])
650 const ScalarT scatterField = (scatterFields_[fieldIndex])(worksetCellIndex,basisId);
653 (*r)[row] = scatterField.val();
654 if(dirichletCounter_!=Teuchos::null) {
656 (*dirichletCounter_)[row] = 1.0;
660 std::vector<double> jacRow(scatterField.size(),0.0);
662 if(!preserveDiagonal_) {
663 int err = Jac->ReplaceMyValues(row, gidCount, scatterField.dx(),
665 TEUCHOS_ASSERT(err==0);
panzer::Traits::Tangent::ScalarT ScalarT
panzer::Traits::Jacobian::ScalarT ScalarT
const Teuchos::RCP< Epetra_Vector > get_f() const
Pushes residual values into the residual vector for a Newton-based solve.