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" 61 #include "Phalanx_DataLayout_MDALayout.hpp" 63 #include "Thyra_SpmdVectorBase.hpp" 64 #include "Thyra_ProductVectorBase.hpp" 65 #include "Thyra_DefaultProductVector.hpp" 66 #include "Thyra_BlockedLinearOpBase.hpp" 67 #include "Thyra_get_Epetra_Operator.hpp" 69 #include "Teuchos_FancyOStream.hpp" 71 #include <unordered_map> 78 template<
typename TRAITS,
typename LO,
typename GO>
82 const Teuchos::ParameterList& p,
83 bool useDiscreteAdjoint)
84 : rowIndexers_(rIndexers)
85 , colIndexers_(cIndexers)
86 , globalDataKey_(
"Residual Scatter Container")
88 std::string scatterName = p.get<std::string>(
"Scatter Name");
90 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
93 const std::vector<std::string>& names =
94 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
97 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
100 scatterIC_ = p.isParameter(
"Scatter Initial Condition") ? p.get<
bool>(
"Scatter Initial Condition") :
false;
102 Teuchos::RCP<PHX::DataLayout> dl = (!scatterIC_) ?
103 p.get< Teuchos::RCP<panzer::PureBasis> >(
"Basis")->functional :
104 p.get< Teuchos::RCP<const panzer::PureBasis> >(
"Basis")->functional;
106 side_subcell_dim_ = p.get<
int>(
"Side Subcell Dimension");
107 local_side_id_ = p.get<
int>(
"Local Side ID");
112 for (std::size_t eq = 0; eq < names.size(); ++eq) {
113 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
119 checkApplyBC_ = p.isParameter(
"Check Apply BC") ? p.get<
bool>(
"Check Apply BC") :
false;
121 applyBC_.resize(names.size());
122 for (std::size_t eq = 0; eq < names.size(); ++eq) {
123 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string(
"APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
124 this->addDependentField(applyBC_[eq]);
129 this->addEvaluatedField(*scatterHolder_);
131 if (p.isType<std::string>(
"Global Data Key"))
132 globalDataKey_ = p.get<std::string>(
"Global Data Key");
134 this->setName(scatterName+
" Scatter Residual");
138 template<
typename TRAITS,
typename LO,
typename GO>
149 std::string fieldName = fieldMap_->find(
scatterFields_[fd].fieldTag().name())->second;
152 subFieldIds_[fd] = rowIndexers_[indexerIds_[fd]]->getFieldNum(fieldName);
158 this->utils.setFieldData(applyBC_[fd],fm);
166 template<
typename TRAITS,
typename LO,
typename GO>
173 using Teuchos::rcp_dynamic_cast;
177 Teuchos::RCP<BLOC> blockContainer
178 = Teuchos::rcp_dynamic_cast<BLOC>(d.gedc.getDataObject(
"Dirichlet Counter"),
true);
181 TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
184 Teuchos::RCP<const BLOC> blockedContainer = Teuchos::rcp_dynamic_cast<
const BLOC>(d.gedc.getDataObject(globalDataKey_));
185 Teuchos::RCP<const ELOC> epetraContainer = Teuchos::rcp_dynamic_cast<
const ELOC>(d.gedc.getDataObject(globalDataKey_));
188 if(blockedContainer!=Teuchos::null)
190 rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer->get_f(),
true) :
191 rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer->get_x(),
true);
192 else if(epetraContainer!=Teuchos::null)
194 Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_f_th()) :
195 Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_x_th());
197 TEUCHOS_ASSERT(r_!=Teuchos::null);
201 template<
typename TRAITS,
typename LO,
typename GO>
206 using Teuchos::ArrayRCP;
207 using Teuchos::ptrFromRef;
208 using Teuchos::rcp_dynamic_cast;
211 using Thyra::SpmdVectorBase;
215 std::string blockId = this->wda(workset).block_id;
216 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
224 Teuchos::ArrayRCP<double> local_r;
225 Teuchos::ArrayRCP<double> local_dc;
226 for(std::size_t fieldIndex = 0; fieldIndex <
scatterFields_.size(); fieldIndex++) {
227 int rowIndexer = indexerIds_[fieldIndex];
228 int subFieldNum = subFieldIds_[fieldIndex];
230 rcp_dynamic_cast<SpmdVectorBase<double> >(dirichletCounter_->getNonconstVectorBlock(rowIndexer))
231 ->getNonconstLocalData(ptrFromRef(local_dc));
234 rcp_dynamic_cast<SpmdVectorBase<double> >(r_->getNonconstVectorBlock(rowIndexer))
235 ->getNonconstLocalData(ptrFromRef(local_r));
237 auto subRowIndexer = rowIndexers_[rowIndexer];
240 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
241 std::size_t cellLocalId = localCellIds[worksetCellIndex];
243 const std::vector<LO> & LIDs = subRowIndexer->getElementLIDs(cellLocalId);
247 const std::pair<std::vector<int>,std::vector<int> > & indicePair
248 = subRowIndexer->getGIDFieldOffsets_closure(blockId,subFieldNum, side_subcell_dim_, local_side_id_);
249 const std::vector<int> & elmtOffset = indicePair.first;
250 const std::vector<int> & basisIdMap = indicePair.second;
259 int basisId = basisIdMap[
basis];
262 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
265 local_r[lid] = (
scatterFields_[fieldIndex])(worksetCellIndex,basisId);
272 const std::vector<int> & elmtOffset = subRowIndexer->getGIDFieldOffsets(blockId,subFieldNum);
296 template<
typename TRAITS,
typename LO,
typename GO>
300 const Teuchos::ParameterList& p,
301 bool useDiscreteAdjoint)
302 : rowIndexers_(rIndexers)
303 , colIndexers_(cIndexers)
304 , globalDataKey_(
"Residual Scatter Container")
306 std::string scatterName = p.get<std::string>(
"Scatter Name");
308 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
311 const std::vector<std::string>& names =
312 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
315 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
318 scatterIC_ = p.isParameter(
"Scatter Initial Condition") ? p.get<
bool>(
"Scatter Initial Condition") :
false;
320 Teuchos::RCP<PHX::DataLayout> dl = (!scatterIC_) ?
321 p.get< Teuchos::RCP<panzer::PureBasis> >(
"Basis")->functional :
322 p.get< Teuchos::RCP<const panzer::PureBasis> >(
"Basis")->functional;
324 side_subcell_dim_ = p.get<
int>(
"Side Subcell Dimension");
325 local_side_id_ = p.get<
int>(
"Local Side ID");
330 for (std::size_t eq = 0; eq < names.size(); ++eq) {
331 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
337 checkApplyBC_ = p.isParameter(
"Check Apply BC") ? p.get<
bool>(
"Check Apply BC") :
false;
339 applyBC_.resize(names.size());
340 for (std::size_t eq = 0; eq < names.size(); ++eq) {
341 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string(
"APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
342 this->addDependentField(applyBC_[eq]);
347 this->addEvaluatedField(*scatterHolder_);
349 if (p.isType<std::string>(
"Global Data Key"))
350 globalDataKey_ = p.get<std::string>(
"Global Data Key");
352 this->setName(scatterName+
" Scatter Tangent");
356 template<
typename TRAITS,
typename LO,
typename GO>
367 std::string fieldName = fieldMap_->find(
scatterFields_[fd].fieldTag().name())->second;
370 subFieldIds_[fd] = rowIndexers_[indexerIds_[fd]]->getFieldNum(fieldName);
376 this->utils.setFieldData(applyBC_[fd],fm);
384 template<
typename TRAITS,
typename LO,
typename GO>
391 using Teuchos::rcp_dynamic_cast;
395 Teuchos::RCP<BLOC> blockContainer
396 = Teuchos::rcp_dynamic_cast<BLOC>(d.gedc.getDataObject(
"Dirichlet Counter"),
true);
399 TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
402 Teuchos::RCP<const BLOC> blockedContainer = Teuchos::rcp_dynamic_cast<
const BLOC>(d.gedc.getDataObject(globalDataKey_));
403 Teuchos::RCP<const ELOC> epetraContainer = Teuchos::rcp_dynamic_cast<
const ELOC>(d.gedc.getDataObject(globalDataKey_));
406 if(blockedContainer!=Teuchos::null)
408 rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer->get_f(),
true) :
409 rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer->get_x(),
true);
410 else if(epetraContainer!=Teuchos::null)
412 Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_f_th()) :
413 Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_x_th());
415 TEUCHOS_ASSERT(r_!=Teuchos::null);
419 template<
typename TRAITS,
typename LO,
typename GO>
423 TEUCHOS_ASSERT(
false);
426 using Teuchos::ArrayRCP;
427 using Teuchos::ptrFromRef;
428 using Teuchos::rcp_dynamic_cast;
431 using Thyra::SpmdVectorBase;
435 std::string blockId = this->wda(workset).block_id;
436 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
444 Teuchos::ArrayRCP<double> local_r;
445 Teuchos::ArrayRCP<double> local_dc;
446 for(std::size_t fieldIndex = 0; fieldIndex <
scatterFields_.size(); fieldIndex++) {
447 int rowIndexer = indexerIds_[fieldIndex];
448 int subFieldNum = subFieldIds_[fieldIndex];
450 rcp_dynamic_cast<SpmdVectorBase<double> >(dirichletCounter_->getNonconstVectorBlock(rowIndexer))
451 ->getNonconstLocalData(ptrFromRef(local_dc));
454 rcp_dynamic_cast<SpmdVectorBase<double> >(r_->getNonconstVectorBlock(rowIndexer))
455 ->getNonconstLocalData(ptrFromRef(local_r));
457 auto subRowIndexer = rowIndexers_[rowIndexer];
460 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
461 std::size_t cellLocalId = localCellIds[worksetCellIndex];
463 const std::vector<LO> & LIDs = subRowIndexer->getElementLIDs(cellLocalId);
467 const std::pair<std::vector<int>,std::vector<int> > & indicePair
468 = subRowIndexer->getGIDFieldOffsets_closure(blockId,subFieldNum, side_subcell_dim_, local_side_id_);
469 const std::vector<int> & elmtOffset = indicePair.first;
470 const std::vector<int> & basisIdMap = indicePair.second;
479 int basisId = basisIdMap[
basis];
482 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
485 local_r[lid] = (
scatterFields_[fieldIndex])(worksetCellIndex,basisId).val();
492 const std::vector<int> & elmtOffset = subRowIndexer->getGIDFieldOffsets(blockId,subFieldNum);
515 template<
typename TRAITS,
typename LO,
typename GO>
519 const Teuchos::ParameterList& p,
520 bool useDiscreteAdjoint)
521 : rowIndexers_(rIndexers)
522 , colIndexers_(cIndexers)
523 , globalDataKey_(
"Residual Scatter Container")
525 std::string scatterName = p.get<std::string>(
"Scatter Name");
527 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
530 const std::vector<std::string>& names =
531 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
534 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
536 Teuchos::RCP<PHX::DataLayout> dl =
537 p.get< Teuchos::RCP<panzer::PureBasis> >(
"Basis")->functional;
539 side_subcell_dim_ = p.get<
int>(
"Side Subcell Dimension");
540 local_side_id_ = p.get<
int>(
"Local Side ID");
544 for (std::size_t eq = 0; eq < names.size(); ++eq) {
545 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
551 checkApplyBC_ = p.get<
bool>(
"Check Apply BC");
553 applyBC_.resize(names.size());
554 for (std::size_t eq = 0; eq < names.size(); ++eq) {
555 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string(
"APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
556 this->addDependentField(applyBC_[eq]);
561 this->addEvaluatedField(*scatterHolder_);
563 if (p.isType<std::string>(
"Global Data Key"))
564 globalDataKey_ = p.get<std::string>(
"Global Data Key");
566 if(colIndexers_.size()==0)
567 colIndexers_ = rowIndexers_;
569 this->setName(scatterName+
" Scatter Residual (Jacobian)");
573 template<
typename TRAITS,
typename LO,
typename GO>
584 std::string fieldName = fieldMap_->find(
scatterFields_[fd].fieldTag().name())->second;
587 subFieldIds_[fd] = rowIndexers_[indexerIds_[fd]]->getFieldNum(fieldName);
593 this->utils.setFieldData(applyBC_[fd],fm);
602 template<
typename TRAITS,
typename LO,
typename GO>
608 using Teuchos::rcp_dynamic_cast;
611 Teuchos::RCP<const BLOC> blockContainer
612 = rcp_dynamic_cast<
const BLOC>(d.gedc.getDataObject(
"Dirichlet Counter"),
true);
615 TEUCHOS_ASSERT(!Teuchos::is_null(dirichletCounter_));
618 blockContainer = rcp_dynamic_cast<
const BLOC>(d.gedc.getDataObject(globalDataKey_),
true);
619 TEUCHOS_ASSERT(!Teuchos::is_null(blockContainer));
626 template<
typename TRAITS,
typename LO,
typename GO>
631 using Teuchos::ArrayRCP;
632 using Teuchos::ptrFromRef;
633 using Teuchos::rcp_dynamic_cast;
635 using Thyra::SpmdVectorBase;
638 std::string blockId = this->wda(workset).block_id;
639 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
641 int numFieldBlocks = Teuchos::as<int>(colIndexers_.size());
643 std::vector<int> blockOffsets;
646 std::unordered_map<std::pair<int,int>,Teuchos::RCP<Epetra_CrsMatrix>,
panzer::pair_hash> jacEpetraBlocks;
653 for(std::size_t fieldIndex = 0; fieldIndex <
scatterFields_.size(); fieldIndex++) {
654 int rowIndexer = indexerIds_[fieldIndex];
655 int subFieldNum = subFieldIds_[fieldIndex];
658 Teuchos::ArrayRCP<double> local_dc;
659 rcp_dynamic_cast<SpmdVectorBase<double> >(dirichletCounter_->getNonconstVectorBlock(rowIndexer))
660 ->getNonconstLocalData(ptrFromRef(local_dc));
663 Teuchos::ArrayRCP<double> local_r;
664 if(r_!=Teuchos::null)
665 rcp_dynamic_cast<SpmdVectorBase<double> >(r_->getNonconstVectorBlock(rowIndexer))
666 ->getNonconstLocalData(ptrFromRef(local_r));
668 auto subRowIndexer = rowIndexers_[rowIndexer];
669 auto subIndicePair = subRowIndexer->getGIDFieldOffsets_closure(blockId,subFieldNum, side_subcell_dim_, local_side_id_);
670 const std::vector<int> & subElmtOffset = subIndicePair.first;
671 const std::vector<int> & subBasisIdMap = subIndicePair.second;
674 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
675 std::size_t cellLocalId = localCellIds[worksetCellIndex];
677 const std::vector<LO> & rLIDs = subRowIndexer->getElementLIDs(cellLocalId);
686 int basisId = subBasisIdMap[
basis];
689 if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
693 for(
int colIndexer=0;colIndexer<numFieldBlocks;colIndexer++) {
694 int start = blockOffsets[colIndexer];
695 int end = blockOffsets[colIndexer+1];
701 std::pair<int,int> blockIndex = std::make_pair(rowIndexer,colIndexer);
702 Teuchos::RCP<Epetra_CrsMatrix> subJac = jacEpetraBlocks[blockIndex];
704 if(subJac==Teuchos::null) {
705 Teuchos::RCP<Thyra::LinearOpBase<double> > tOp = Jac_->getNonconstBlock(blockIndex.first,blockIndex.second);
708 if(Teuchos::is_null(tOp))
711 Teuchos::RCP<Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*tOp);
712 subJac = rcp_dynamic_cast<Epetra_CrsMatrix>(eOp,
true);
713 jacEpetraBlocks[blockIndex] = subJac;
717 int * rowIndices = 0;
718 double * rowValues = 0;
720 subJac->ExtractMyRowView(lid,numEntries,rowValues,rowIndices);
722 for(
int i=0;i<numEntries;i++)
728 if(r_!=Teuchos::null)
729 local_r[lid] = scatterField.val();
733 std::vector<double> jacRow(scatterField.size(),0.0);
735 for(
int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
736 jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex);
738 for(
int colIndexer=0;colIndexer<numFieldBlocks;colIndexer++) {
739 int start = blockOffsets[colIndexer];
740 int end = blockOffsets[colIndexer+1];
745 auto subColIndexer = colIndexers_[colIndexer];
746 const std::vector<LO> & cLIDs = subColIndexer->getElementLIDs(cellLocalId);
748 TEUCHOS_ASSERT(end-start==Teuchos::as<int>(cLIDs.size()));
751 std::pair<int,int> blockIndex = std::make_pair(rowIndexer,colIndexer);
752 Teuchos::RCP<Epetra_CrsMatrix> subJac = jacEpetraBlocks[blockIndex];
755 if(subJac==Teuchos::null) {
756 Teuchos::RCP<Thyra::LinearOpBase<double> > tOp = Jac_->getNonconstBlock(blockIndex.first,blockIndex.second);
759 if(Teuchos::is_null(tOp))
762 Teuchos::RCP<Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*tOp);
763 subJac = rcp_dynamic_cast<Epetra_CrsMatrix>(eOp,
true);
764 jacEpetraBlocks[blockIndex] = subJac;
768 int err = subJac->ReplaceMyValues(lid, end-start, &jacRow[start],&cLIDs[0]);
770 std::stringstream ss;
771 ss <<
"Failed inserting row: " <<
" (" << lid <<
"): ";
772 for(
int i=0;i<end-start;i++)
773 ss << cLIDs[i] <<
" ";
775 ss <<
"Into block " << rowIndexer <<
", " << colIndexer << std::endl;
777 ss <<
"scatter field = ";
781 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)
std::vector< PHX::MDField< const ScalarT, panzer::Cell, panzer::Point > > scatterFields_
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)
PHX::MDField< ScalarT > vector
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
Teuchos::RCP< const panzer::PureBasis > basis
Interpolates basis DOF values to IP DOF values.