43 #ifndef PANZER_SCATTER_DIRICHLET_RESIDUAL_BLOCEDEPETRA_IMPL_HPP 44 #define PANZER_SCATTER_DIRICHLET_RESIDUAL_BLOCEDEPETRA_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 "Thyra_SpmdVectorBase.hpp" 65 #include "Thyra_ProductVectorBase.hpp" 66 #include "Thyra_DefaultProductVector.hpp" 67 #include "Thyra_BlockedLinearOpBase.hpp" 68 #include "Thyra_get_Epetra_Operator.hpp" 70 #include "Teuchos_FancyOStream.hpp" 72 #include <unordered_map> 79 template<
typename TRAITS,
typename LO,
typename GO>
83 const Teuchos::ParameterList& p,
85 : rowIndexers_(rIndexers)
86 , colIndexers_(cIndexers)
87 , globalDataKey_(
"Residual Scatter Container")
89 std::string scatterName = p.get<std::string>(
"Scatter Name");
91 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
94 const std::vector<std::string>& names =
95 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
98 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
101 scatterIC_ = p.isParameter(
"Scatter Initial Condition") ? p.get<
bool>(
"Scatter Initial Condition") :
false;
103 Teuchos::RCP<PHX::DataLayout> dl = (!scatterIC_) ?
104 p.get< Teuchos::RCP<panzer::PureBasis> >(
"Basis")->functional :
105 p.get< Teuchos::RCP<const panzer::PureBasis> >(
"Basis")->functional;
107 side_subcell_dim_ = p.get<
int>(
"Side Subcell Dimension");
108 local_side_id_ = p.get<
int>(
"Local Side ID");
112 scatterFields_.resize(names.size());
113 for (std::size_t eq = 0; eq < names.size(); ++eq) {
114 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
117 this->addDependentField(scatterFields_[eq]);
120 checkApplyBC_ = p.isParameter(
"Check Apply BC") ? p.get<
bool>(
"Check Apply BC") :
false;
122 applyBC_.resize(names.size());
123 for (std::size_t eq = 0; eq < names.size(); ++eq) {
124 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string(
"APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
125 this->addDependentField(applyBC_[eq]);
130 this->addEvaluatedField(*scatterHolder_);
132 if (p.isType<std::string>(
"Global Data Key"))
133 globalDataKey_ = p.get<std::string>(
"Global Data Key");
135 this->setName(scatterName+
" Scatter Residual");
139 template<
typename TRAITS,
typename LO,
typename GO>
144 indexerIds_.resize(scatterFields_.size());
145 subFieldIds_.resize(scatterFields_.size());
148 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
150 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
153 subFieldIds_[fd] = rowIndexers_[indexerIds_[fd]]->getFieldNum(fieldName);
157 num_nodes = scatterFields_[0].extent(1);
161 template<
typename TRAITS,
typename LO,
typename GO>
168 using Teuchos::rcp_dynamic_cast;
172 Teuchos::RCP<BLOC> blockContainer
173 = Teuchos::rcp_dynamic_cast<BLOC>(d.gedc->getDataObject(
"Dirichlet Counter"),
true);
176 TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
179 Teuchos::RCP<const BLOC> blockedContainer = Teuchos::rcp_dynamic_cast<
const BLOC>(d.gedc->getDataObject(globalDataKey_));
180 Teuchos::RCP<const ELOC> epetraContainer = Teuchos::rcp_dynamic_cast<
const ELOC>(d.gedc->getDataObject(globalDataKey_));
183 if(blockedContainer!=Teuchos::null)
185 rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer->get_f(),
true) :
186 rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer->get_x(),
true);
187 else if(epetraContainer!=Teuchos::null)
189 Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_f_th()) :
190 Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_x_th());
192 TEUCHOS_ASSERT(r_!=Teuchos::null);
196 template<
typename TRAITS,
typename LO,
typename GO>
201 using Teuchos::ArrayRCP;
202 using Teuchos::ptrFromRef;
203 using Teuchos::rcp_dynamic_cast;
206 using Thyra::SpmdVectorBase;
210 std::string blockId = this->wda(workset).block_id;
211 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
219 Teuchos::ArrayRCP<double> local_r;
220 Teuchos::ArrayRCP<double> local_dc;
221 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
222 int rowIndexer = indexerIds_[fieldIndex];
223 int subFieldNum = subFieldIds_[fieldIndex];
225 rcp_dynamic_cast<SpmdVectorBase<double> >(dirichletCounter_->getNonconstVectorBlock(rowIndexer))
226 ->getNonconstLocalData(ptrFromRef(local_dc));
229 rcp_dynamic_cast<SpmdVectorBase<double> >(r_->getNonconstVectorBlock(rowIndexer))
230 ->getNonconstLocalData(ptrFromRef(local_r));
232 auto subRowIndexer = rowIndexers_[rowIndexer];
235 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
236 std::size_t cellLocalId = localCellIds[worksetCellIndex];
238 auto LIDs = subRowIndexer->getElementLIDs(cellLocalId);
242 const std::pair<std::vector<int>,std::vector<int> > & indicePair
243 = subRowIndexer->getGIDFieldOffsets_closure(blockId,subFieldNum, side_subcell_dim_, local_side_id_);
244 const std::vector<int> & elmtOffset = indicePair.first;
245 const std::vector<int> & basisIdMap = indicePair.second;
248 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
249 int offset = elmtOffset[basis];
250 int lid = LIDs[offset];
254 int basisId = basisIdMap[basis];
257 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
260 local_r[lid] = (scatterFields_[fieldIndex])(worksetCellIndex,basisId);
267 const std::vector<int> & elmtOffset = subRowIndexer->getGIDFieldOffsets(blockId,subFieldNum);
270 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
271 int offset = elmtOffset[basis];
272 int lid = LIDs[offset];
276 local_r[lid] = (scatterFields_[fieldIndex])(worksetCellIndex,basis);
291 template<
typename TRAITS,
typename LO,
typename GO>
295 const Teuchos::ParameterList& p,
297 : rowIndexers_(rIndexers)
298 , colIndexers_(cIndexers)
299 , globalDataKey_(
"Residual Scatter Container")
301 std::string scatterName = p.get<std::string>(
"Scatter Name");
303 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
306 const std::vector<std::string>& names =
307 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
310 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
313 scatterIC_ = p.isParameter(
"Scatter Initial Condition") ? p.get<
bool>(
"Scatter Initial Condition") :
false;
315 Teuchos::RCP<PHX::DataLayout> dl = (!scatterIC_) ?
316 p.get< Teuchos::RCP<panzer::PureBasis> >(
"Basis")->functional :
317 p.get< Teuchos::RCP<const panzer::PureBasis> >(
"Basis")->functional;
319 side_subcell_dim_ = p.get<
int>(
"Side Subcell Dimension");
320 local_side_id_ = p.get<
int>(
"Local Side ID");
324 scatterFields_.resize(names.size());
325 for (std::size_t eq = 0; eq < names.size(); ++eq) {
326 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
329 this->addDependentField(scatterFields_[eq]);
332 checkApplyBC_ = p.isParameter(
"Check Apply BC") ? p.get<
bool>(
"Check Apply BC") :
false;
334 applyBC_.resize(names.size());
335 for (std::size_t eq = 0; eq < names.size(); ++eq) {
336 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string(
"APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
337 this->addDependentField(applyBC_[eq]);
342 this->addEvaluatedField(*scatterHolder_);
344 if (p.isType<std::string>(
"Global Data Key"))
345 globalDataKey_ = p.get<std::string>(
"Global Data Key");
347 this->setName(scatterName+
" Scatter Tangent");
351 template<
typename TRAITS,
typename LO,
typename GO>
356 indexerIds_.resize(scatterFields_.size());
357 subFieldIds_.resize(scatterFields_.size());
360 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
362 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
365 subFieldIds_[fd] = rowIndexers_[indexerIds_[fd]]->getFieldNum(fieldName);
369 num_nodes = scatterFields_[0].extent(1);
373 template<
typename TRAITS,
typename LO,
typename GO>
380 using Teuchos::rcp_dynamic_cast;
384 Teuchos::RCP<BLOC> blockContainer
385 = Teuchos::rcp_dynamic_cast<BLOC>(d.gedc->getDataObject(
"Dirichlet Counter"),
true);
388 TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
391 Teuchos::RCP<const BLOC> blockedContainer = Teuchos::rcp_dynamic_cast<
const BLOC>(d.gedc->getDataObject(globalDataKey_));
392 Teuchos::RCP<const ELOC> epetraContainer = Teuchos::rcp_dynamic_cast<
const ELOC>(d.gedc->getDataObject(globalDataKey_));
395 if(blockedContainer!=Teuchos::null)
397 rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer->get_f(),
true) :
398 rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer->get_x(),
true);
399 else if(epetraContainer!=Teuchos::null)
401 Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_f_th()) :
402 Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_x_th());
404 TEUCHOS_ASSERT(r_!=Teuchos::null);
408 template<
typename TRAITS,
typename LO,
typename GO>
412 TEUCHOS_ASSERT(
false);
415 using Teuchos::ArrayRCP;
416 using Teuchos::ptrFromRef;
417 using Teuchos::rcp_dynamic_cast;
420 using Thyra::SpmdVectorBase;
424 std::string blockId = this->wda(workset).block_id;
425 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
433 Teuchos::ArrayRCP<double> local_r;
434 Teuchos::ArrayRCP<double> local_dc;
435 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
436 int rowIndexer = indexerIds_[fieldIndex];
437 int subFieldNum = subFieldIds_[fieldIndex];
439 rcp_dynamic_cast<SpmdVectorBase<double> >(dirichletCounter_->getNonconstVectorBlock(rowIndexer))
440 ->getNonconstLocalData(ptrFromRef(local_dc));
443 rcp_dynamic_cast<SpmdVectorBase<double> >(r_->getNonconstVectorBlock(rowIndexer))
444 ->getNonconstLocalData(ptrFromRef(local_r));
446 auto subRowIndexer = rowIndexers_[rowIndexer];
449 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
450 std::size_t cellLocalId = localCellIds[worksetCellIndex];
452 auto LIDs = subRowIndexer->getElementLIDs(cellLocalId);
456 const std::pair<std::vector<int>,std::vector<int> > & indicePair
457 = subRowIndexer->getGIDFieldOffsets_closure(blockId,subFieldNum, side_subcell_dim_, local_side_id_);
458 const std::vector<int> & elmtOffset = indicePair.first;
459 const std::vector<int> & basisIdMap = indicePair.second;
462 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
463 int offset = elmtOffset[basis];
464 int lid = LIDs[offset];
468 int basisId = basisIdMap[basis];
471 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
474 local_r[lid] = (scatterFields_[fieldIndex])(worksetCellIndex,basisId).val();
481 const std::vector<int> & elmtOffset = subRowIndexer->getGIDFieldOffsets(blockId,subFieldNum);
484 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
485 int offset = elmtOffset[basis];
486 int lid = LIDs[offset];
490 local_r[lid] = (scatterFields_[fieldIndex])(worksetCellIndex,basis).val();
504 template<
typename TRAITS,
typename LO,
typename GO>
508 const Teuchos::ParameterList& p,
510 : rowIndexers_(rIndexers)
511 , colIndexers_(cIndexers)
512 , globalDataKey_(
"Residual Scatter Container")
514 std::string scatterName = p.get<std::string>(
"Scatter Name");
516 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
519 const std::vector<std::string>& names =
520 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
523 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
525 Teuchos::RCP<PHX::DataLayout> dl =
526 p.get< Teuchos::RCP<panzer::PureBasis> >(
"Basis")->functional;
528 side_subcell_dim_ = p.get<
int>(
"Side Subcell Dimension");
529 local_side_id_ = p.get<
int>(
"Local Side ID");
532 scatterFields_.resize(names.size());
533 for (std::size_t eq = 0; eq < names.size(); ++eq) {
534 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
537 this->addDependentField(scatterFields_[eq]);
540 checkApplyBC_ = p.get<
bool>(
"Check Apply BC");
542 applyBC_.resize(names.size());
543 for (std::size_t eq = 0; eq < names.size(); ++eq) {
544 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string(
"APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
545 this->addDependentField(applyBC_[eq]);
550 this->addEvaluatedField(*scatterHolder_);
552 if (p.isType<std::string>(
"Global Data Key"))
553 globalDataKey_ = p.get<std::string>(
"Global Data Key");
555 if(colIndexers_.size()==0)
556 colIndexers_ = rowIndexers_;
558 this->setName(scatterName+
" Scatter Residual (Jacobian)");
562 template<
typename TRAITS,
typename LO,
typename GO>
567 indexerIds_.resize(scatterFields_.size());
568 subFieldIds_.resize(scatterFields_.size());
571 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
573 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
576 subFieldIds_[fd] = rowIndexers_[indexerIds_[fd]]->getFieldNum(fieldName);
580 num_nodes = scatterFields_[0].extent(1);
581 num_eq = scatterFields_.size();
585 template<
typename TRAITS,
typename LO,
typename GO>
591 using Teuchos::rcp_dynamic_cast;
594 Teuchos::RCP<const BLOC> blockContainer
595 = rcp_dynamic_cast<
const BLOC>(d.gedc->getDataObject(
"Dirichlet Counter"),
true);
598 TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
601 blockContainer = rcp_dynamic_cast<
const BLOC>(d.gedc->getDataObject(globalDataKey_),
true);
602 TEUCHOS_ASSERT(!Teuchos::is_null(blockContainer));
609 template<
typename TRAITS,
typename LO,
typename GO>
614 using Teuchos::ArrayRCP;
615 using Teuchos::ptrFromRef;
616 using Teuchos::rcp_dynamic_cast;
618 using Thyra::SpmdVectorBase;
621 std::string blockId = this->wda(workset).block_id;
622 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
624 int numFieldBlocks = Teuchos::as<int>(colIndexers_.size());
626 std::vector<int> blockOffsets;
629 std::unordered_map<std::pair<int,int>,Teuchos::RCP<Epetra_CrsMatrix>,
panzer::pair_hash> jacEpetraBlocks;
636 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
637 int rowIndexer = indexerIds_[fieldIndex];
638 int subFieldNum = subFieldIds_[fieldIndex];
641 Teuchos::ArrayRCP<double> local_dc;
642 rcp_dynamic_cast<SpmdVectorBase<double> >(dirichletCounter_->getNonconstVectorBlock(rowIndexer))
643 ->getNonconstLocalData(ptrFromRef(local_dc));
646 Teuchos::ArrayRCP<double> local_r;
647 if(r_!=Teuchos::null)
648 rcp_dynamic_cast<SpmdVectorBase<double> >(r_->getNonconstVectorBlock(rowIndexer))
649 ->getNonconstLocalData(ptrFromRef(local_r));
651 auto subRowIndexer = rowIndexers_[rowIndexer];
652 auto subIndicePair = subRowIndexer->getGIDFieldOffsets_closure(blockId,subFieldNum, side_subcell_dim_, local_side_id_);
653 const std::vector<int> & subElmtOffset = subIndicePair.first;
654 const std::vector<int> & subBasisIdMap = subIndicePair.second;
657 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
658 std::size_t cellLocalId = localCellIds[worksetCellIndex];
660 auto rLIDs = subRowIndexer->getElementLIDs(cellLocalId);
663 for(std::size_t basis=0;basis<subElmtOffset.size();basis++) {
664 int offset = subElmtOffset[basis];
665 int lid = rLIDs[offset];
669 int basisId = subBasisIdMap[basis];
672 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
676 for(
int colIndexer=0;colIndexer<numFieldBlocks;colIndexer++) {
677 int start = blockOffsets[colIndexer];
678 int end = blockOffsets[colIndexer+1];
684 std::pair<int,int> blockIndex = std::make_pair(rowIndexer,colIndexer);
685 Teuchos::RCP<Epetra_CrsMatrix> subJac = jacEpetraBlocks[blockIndex];
687 if(subJac==Teuchos::null) {
688 Teuchos::RCP<Thyra::LinearOpBase<double> > tOp = Jac_->getNonconstBlock(blockIndex.first,blockIndex.second);
691 if(Teuchos::is_null(tOp))
694 Teuchos::RCP<Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*tOp);
695 subJac = rcp_dynamic_cast<Epetra_CrsMatrix>(eOp,
true);
696 jacEpetraBlocks[blockIndex] = subJac;
700 int * rowIndices = 0;
701 double * rowValues = 0;
703 subJac->ExtractMyRowView(lid,numEntries,rowValues,rowIndices);
705 for(
int i=0;i<numEntries;i++)
709 const ScalarT scatterField = (scatterFields_[fieldIndex])(worksetCellIndex,basisId);
711 if(r_!=Teuchos::null)
712 local_r[lid] = scatterField.val();
716 std::vector<double> jacRow(scatterField.size(),0.0);
718 for(
int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
719 jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex);
721 for(
int colIndexer=0;colIndexer<numFieldBlocks;colIndexer++) {
722 int start = blockOffsets[colIndexer];
723 int end = blockOffsets[colIndexer+1];
728 auto subColIndexer = colIndexers_[colIndexer];
729 auto cLIDs = subColIndexer->getElementLIDs(cellLocalId);
731 TEUCHOS_ASSERT(end-start==Teuchos::as<int>(cLIDs.size()));
734 std::pair<int,int> blockIndex = std::make_pair(rowIndexer,colIndexer);
735 Teuchos::RCP<Epetra_CrsMatrix> subJac = jacEpetraBlocks[blockIndex];
738 if(subJac==Teuchos::null) {
739 Teuchos::RCP<Thyra::LinearOpBase<double> > tOp = Jac_->getNonconstBlock(blockIndex.first,blockIndex.second);
742 if(Teuchos::is_null(tOp))
745 Teuchos::RCP<Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*tOp);
746 subJac = rcp_dynamic_cast<Epetra_CrsMatrix>(eOp,
true);
747 jacEpetraBlocks[blockIndex] = subJac;
751 int err = subJac->ReplaceMyValues(lid, end-start, &jacRow[start],&cLIDs[0]);
753 std::stringstream ss;
754 ss <<
"Failed inserting row: " <<
" (" << lid <<
"): ";
755 for(
int i=0;i<end-start;i++)
756 ss << cLIDs[i] <<
" ";
758 ss <<
"Into block " << rowIndexer <<
", " << colIndexer << std::endl;
760 ss <<
"scatter field = ";
761 scatterFields_[fieldIndex].print(ss);
764 TEUCHOS_TEST_FOR_EXCEPTION(err!=0,std::runtime_error,ss.str());
int getFieldBlock(const std::string &fieldName, const std::vector< Teuchos::RCP< UniqueGlobalIndexer< LocalOrdinalT, GlobalOrdinalT > > > &ugis)
void computeBlockOffsets(const std::string &blockId, const std::vector< Teuchos::RCP< UniqueGlobalIndexer< LocalOrdinalT, GlobalOrdinalT > > > &ugis, std::vector< int > &blockOffsets)
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &vm)
Pushes residual values into the residual vector for a Newton-based solve.
void evaluateFields(typename TRAITS::EvalData d)
ScatterDirichletResidual_BlockedEpetra(const Teuchos::ParameterList &p)
panzer::Traits::Jacobian::ScalarT ScalarT