43 #ifndef PANZER_SCATTER_DIRICHLET_RESIDUAL_TPETRA_IMPL_HPP 44 #define PANZER_SCATTER_DIRICHLET_RESIDUAL_TPETRA_IMPL_HPP 46 #include "Teuchos_RCP.hpp" 47 #include "Teuchos_Assert.hpp" 49 #include "Phalanx_DataLayout.hpp" 58 #include "Phalanx_DataLayout_MDALayout.hpp" 60 #include "Teuchos_FancyOStream.hpp" 67 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
70 const Teuchos::ParameterList& p)
71 : globalIndexer_(indexer)
72 , globalDataKey_(
"Residual Scatter Container")
74 std::string scatterName = p.get<std::string>(
"Scatter Name");
76 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
79 const std::vector<std::string>& names =
80 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
83 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
86 scatterIC_ = p.isParameter(
"Scatter Initial Condition") ? p.get<
bool>(
"Scatter Initial Condition") :
false;
88 Teuchos::RCP<PHX::DataLayout> dl = (!scatterIC_) ?
89 p.get< Teuchos::RCP<panzer::PureBasis> >(
"Basis")->functional :
90 p.get< Teuchos::RCP<const panzer::PureBasis> >(
"Basis")->functional ;
92 side_subcell_dim_ = p.get<
int>(
"Side Subcell Dimension");
93 local_side_id_ = p.get<
int>(
"Local Side ID");
97 scatterFields_.resize(names.size());
98 for (std::size_t eq = 0; eq < names.size(); ++eq) {
99 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
102 this->addDependentField(scatterFields_[eq]);
105 checkApplyBC_ = p.isParameter(
"Check Apply BC") ? p.get<
bool>(
"Check Apply BC") :
false;
107 applyBC_.resize(names.size());
108 for (std::size_t eq = 0; eq < names.size(); ++eq) {
109 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string(
"APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
110 this->addDependentField(applyBC_[eq]);
115 this->addEvaluatedField(*scatterHolder_);
117 if (p.isType<std::string>(
"Global Data Key"))
118 globalDataKey_ = p.get<std::string>(
"Global Data Key");
120 this->setName(scatterName+
" Scatter Residual");
124 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
129 fieldIds_.resize(scatterFields_.size());
132 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
134 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
135 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
139 num_nodes = scatterFields_[0].extent(1);
143 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
148 tpetraContainer_ = Teuchos::rcp_dynamic_cast<
LOC>(d.gedc->getDataObject(globalDataKey_));
150 if(tpetraContainer_==Teuchos::null) {
152 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<
LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),
true)->getGhostedLOC();
153 tpetraContainer_ = Teuchos::rcp_dynamic_cast<
LOC>(loc);
155 dirichletCounter_ = Teuchos::null;
159 Teuchos::RCP<LOC> tpetraContainer
160 = Teuchos::rcp_dynamic_cast<
LOC>(d.gedc->getDataObject(
"Dirichlet Counter"),
true);
162 dirichletCounter_ = tpetraContainer->
get_f();
163 TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
168 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
172 std::vector<GO> GIDs;
173 std::vector<LO> LIDs;
176 std::string blockId = this->wda(workset).block_id;
177 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
179 Teuchos::RCP<typename LOC::VectorType> r = (!scatterIC_) ?
180 tpetraContainer_->get_f() :
181 tpetraContainer_->get_x();
183 Teuchos::ArrayRCP<double> r_array = r->get1dViewNonConst();
184 Teuchos::ArrayRCP<double> dc_array = dirichletCounter_->get1dViewNonConst();
193 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
194 std::size_t cellLocalId = localCellIds[worksetCellIndex];
196 globalIndexer_->getElementGIDs(cellLocalId,GIDs);
199 LIDs.resize(GIDs.size());
200 for(std::size_t i=0;i<GIDs.size();i++)
201 LIDs[i] = r->getMap()->getLocalElement(GIDs[i]);
204 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
205 int fieldNum = fieldIds_[fieldIndex];
209 const std::pair<std::vector<int>,std::vector<int> > & indicePair
210 = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldNum, side_subcell_dim_, local_side_id_);
211 const std::vector<int> & elmtOffset = indicePair.first;
212 const std::vector<int> & basisIdMap = indicePair.second;
215 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
216 int offset = elmtOffset[basis];
217 LO lid = LIDs[offset];
221 int basisId = basisIdMap[basis];
224 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
227 r_array[lid] = (scatterFields_[fieldIndex])(worksetCellIndex,basisId);
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 LO lid = LIDs[offset];
243 r_array[lid] = (scatterFields_[fieldIndex])(worksetCellIndex,basis);
258 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
261 const Teuchos::ParameterList& p)
262 : globalIndexer_(indexer)
263 , globalDataKey_(
"Residual Scatter Container")
265 std::string scatterName = p.get<std::string>(
"Scatter Name");
267 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
270 const std::vector<std::string>& names =
271 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
274 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
277 scatterIC_ = p.isParameter(
"Scatter Initial Condition") ? p.get<
bool>(
"Scatter Initial Condition") :
false;
279 Teuchos::RCP<PHX::DataLayout> dl = (!scatterIC_) ?
280 p.get< Teuchos::RCP<panzer::PureBasis> >(
"Basis")->functional :
281 p.get< Teuchos::RCP<const panzer::PureBasis> >(
"Basis")->functional ;
283 side_subcell_dim_ = p.get<
int>(
"Side Subcell Dimension");
284 local_side_id_ = p.get<
int>(
"Local Side ID");
288 scatterFields_.resize(names.size());
289 for (std::size_t eq = 0; eq < names.size(); ++eq) {
290 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
293 this->addDependentField(scatterFields_[eq]);
296 checkApplyBC_ = p.isParameter(
"Check Apply BC") ? p.get<
bool>(
"Check Apply BC") :
false;
298 applyBC_.resize(names.size());
299 for (std::size_t eq = 0; eq < names.size(); ++eq) {
300 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string(
"APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
301 this->addDependentField(applyBC_[eq]);
306 this->addEvaluatedField(*scatterHolder_);
308 if (p.isType<std::string>(
"Global Data Key"))
309 globalDataKey_ = p.get<std::string>(
"Global Data Key");
311 this->setName(scatterName+
" Scatter Tangent");
315 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
320 fieldIds_.resize(scatterFields_.size());
323 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
325 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
326 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
330 num_nodes = scatterFields_[0].extent(1);
334 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
339 tpetraContainer_ = Teuchos::rcp_dynamic_cast<
LOC>(d.gedc->getDataObject(globalDataKey_));
341 if(tpetraContainer_==Teuchos::null) {
343 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<
LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),
true)->getGhostedLOC();
344 tpetraContainer_ = Teuchos::rcp_dynamic_cast<
LOC>(loc);
346 dirichletCounter_ = Teuchos::null;
350 Teuchos::RCP<LOC> tpetraContainer
351 = Teuchos::rcp_dynamic_cast<
LOC>(d.gedc->getDataObject(
"Dirichlet Counter"),
true);
353 dirichletCounter_ = tpetraContainer->
get_f();
354 TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
358 using Teuchos::rcp_dynamic_cast;
361 std::vector<std::string> activeParameters =
365 TEUCHOS_ASSERT(!scatterIC_);
366 dfdp_vectors_.clear();
367 for(std::size_t i=0;i<activeParameters.size();i++) {
368 RCP<typename LOC::VectorType> vec =
369 rcp_dynamic_cast<
LOC>(d.gedc->getDataObject(activeParameters[i]),
true)->get_f();
370 Teuchos::ArrayRCP<double> vec_array = vec->get1dViewNonConst();
371 dfdp_vectors_.push_back(vec_array);
376 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
380 std::vector<GO> GIDs;
381 std::vector<LO> LIDs;
384 std::string blockId = this->wda(workset).block_id;
385 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
387 Teuchos::RCP<typename LOC::VectorType> r = (!scatterIC_) ?
388 tpetraContainer_->get_f() :
389 tpetraContainer_->get_x();
391 Teuchos::ArrayRCP<double> r_array = r->get1dViewNonConst();
392 Teuchos::ArrayRCP<double> dc_array = dirichletCounter_->get1dViewNonConst();
401 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
402 std::size_t cellLocalId = localCellIds[worksetCellIndex];
404 globalIndexer_->getElementGIDs(cellLocalId,GIDs);
407 LIDs.resize(GIDs.size());
408 for(std::size_t i=0;i<GIDs.size();i++)
409 LIDs[i] = r->getMap()->getLocalElement(GIDs[i]);
412 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
413 int fieldNum = fieldIds_[fieldIndex];
417 const std::pair<std::vector<int>,std::vector<int> > & indicePair
418 = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldNum, side_subcell_dim_, local_side_id_);
419 const std::vector<int> & elmtOffset = indicePair.first;
420 const std::vector<int> & basisIdMap = indicePair.second;
423 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
424 int offset = elmtOffset[basis];
425 LO lid = LIDs[offset];
429 int basisId = basisIdMap[basis];
432 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
435 ScalarT value = (scatterFields_[fieldIndex])(worksetCellIndex,basisId);
440 for(std::size_t d=0;d<dfdp_vectors_.size();d++)
441 dfdp_vectors_[d][lid] = 0.0;
443 for(
int d=0;d<value.size();d++) {
444 dfdp_vectors_[d][lid] = value.fastAccessDx(d);
452 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
455 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
456 int offset = elmtOffset[basis];
457 LO lid = LIDs[offset];
461 ScalarT value = (scatterFields_[fieldIndex])(worksetCellIndex,basis);
466 for(std::size_t d=0;d<dfdp_vectors_.size();d++)
467 dfdp_vectors_[d][lid] = 0.0;
469 for(
int d=0;d<value.size();d++) {
470 dfdp_vectors_[d][lid] = value.fastAccessDx(d);
485 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
488 const Teuchos::ParameterList& p)
489 : globalIndexer_(indexer)
490 , globalDataKey_(
"Residual Scatter Container")
492 std::string scatterName = p.get<std::string>(
"Scatter Name");
494 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
497 const std::vector<std::string>& names =
498 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
501 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
503 Teuchos::RCP<PHX::DataLayout> dl =
504 p.get< Teuchos::RCP<panzer::PureBasis> >(
"Basis")->functional;
506 side_subcell_dim_ = p.get<
int>(
"Side Subcell Dimension");
507 local_side_id_ = p.get<
int>(
"Local Side ID");
510 scatterFields_.resize(names.size());
511 for (std::size_t eq = 0; eq < names.size(); ++eq) {
512 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
515 this->addDependentField(scatterFields_[eq]);
518 checkApplyBC_ = p.get<
bool>(
"Check Apply BC");
520 applyBC_.resize(names.size());
521 for (std::size_t eq = 0; eq < names.size(); ++eq) {
522 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string(
"APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
523 this->addDependentField(applyBC_[eq]);
528 this->addEvaluatedField(*scatterHolder_);
530 if (p.isType<std::string>(
"Global Data Key"))
531 globalDataKey_ = p.get<std::string>(
"Global Data Key");
533 this->setName(scatterName+
" Scatter Residual (Jacobian)");
537 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
542 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,
typename NodeT>
561 tpetraContainer_ = Teuchos::rcp_dynamic_cast<
LOC>(d.gedc->getDataObject(globalDataKey_));
563 if(tpetraContainer_==Teuchos::null) {
565 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<
LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),
true)->getGhostedLOC();
566 tpetraContainer_ = Teuchos::rcp_dynamic_cast<
LOC>(loc);
568 dirichletCounter_ = Teuchos::null;
572 Teuchos::RCP<LOC> tpetraContainer
573 = Teuchos::rcp_dynamic_cast<
LOC>(d.gedc->getDataObject(
"Dirichlet Counter"),
true);
575 dirichletCounter_ = tpetraContainer->
get_f();
576 TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
581 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
585 std::vector<GO> GIDs;
588 std::string blockId = this->wda(workset).block_id;
589 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
591 Teuchos::RCP<typename LOC::VectorType> r = tpetraContainer_->get_f();
592 Teuchos::RCP<typename LOC::CrsMatrixType> Jac = tpetraContainer_->get_A();
594 Teuchos::ArrayRCP<double> r_array = r->get1dViewNonConst();
595 Teuchos::ArrayRCP<double> dc_array = dirichletCounter_->get1dViewNonConst();
603 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
604 std::size_t cellLocalId = localCellIds[worksetCellIndex];
606 globalIndexer_->getElementGIDs(cellLocalId,GIDs);
607 auto LIDs = globalIndexer_->getElementLIDs(cellLocalId);
610 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
611 int fieldNum = fieldIds_[fieldIndex];
614 const std::pair<std::vector<int>,std::vector<int> > & indicePair
615 = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldNum, side_subcell_dim_, local_side_id_);
616 const std::vector<int> & elmtOffset = indicePair.first;
617 const std::vector<int> & basisIdMap = indicePair.second;
620 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
621 int offset = elmtOffset[basis];
622 LO lid = LIDs[offset];
626 int basisId = basisIdMap[basis];
629 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
634 std::size_t sz = Jac->getNumEntriesInLocalRow(lid);
635 std::size_t numEntries = 0;
636 Teuchos::Array<LO> rowIndices(sz);
637 Teuchos::Array<double> rowValues(sz);
640 Jac->getLocalRowCopy(lid,rowIndices,rowValues,numEntries);
642 for(std::size_t i=0;i<numEntries;i++)
645 Jac->replaceLocalValues(lid,rowIndices,rowValues);
648 GO gid = GIDs[offset];
649 const ScalarT scatterField = (scatterFields_[fieldIndex])(worksetCellIndex,basisId);
651 r_array[lid] = scatterField.val();
655 std::vector<double> jacRow(scatterField.size(),0.0);
657 for(
int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
658 jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex);
659 TEUCHOS_ASSERT(jacRow.size()==GIDs.size());
661 Jac->replaceGlobalValues(gid, GIDs, jacRow);
panzer::Traits::Tangent::ScalarT ScalarT
Pushes residual values into the residual vector for a Newton-based solve.
const Teuchos::RCP< VectorType > get_f() const
panzer::Traits::Jacobian::ScalarT ScalarT