43 #ifndef PANZER_SCATTER_RESIDUAL_BLOCKEDEPETRA_IMPL_HPP 44 #define PANZER_SCATTER_RESIDUAL_BLOCKEDEPETRA_IMPL_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 "Thyra_SpmdVectorBase.hpp" 63 #include "Thyra_ProductVectorBase.hpp" 64 #include "Thyra_DefaultProductVector.hpp" 65 #include "Thyra_BlockedLinearOpBase.hpp" 66 #include "Thyra_DefaultBlockedLinearOp.hpp" 67 #include "Thyra_get_Epetra_Operator.hpp" 69 #include "Phalanx_DataLayout_MDALayout.hpp" 71 #include "Teuchos_FancyOStream.hpp" 77 template<
typename TRAITS,
typename LO,
typename GO>
82 bool useDiscreteAdjoint)
83 : rowIndexers_(rIndexers)
84 , colIndexers_(cIndexers)
85 , globalDataKey_(
"Residual Scatter Container")
87 std::string scatterName = p.
get<std::string>(
"Scatter Name");
92 const std::vector<std::string>& names =
103 for (std::size_t eq = 0; eq < names.size(); ++eq) {
104 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
111 this->addEvaluatedField(*scatterHolder_);
113 if (p.
isType<std::string>(
"Global Data Key"))
114 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
115 if (p.
isType<
bool>(
"Use Discrete Adjoint"))
116 useDiscreteAdjoint = p.
get<
bool>(
"Use Discrete Adjoint");
118 this->setName(scatterName+
" Scatter Residual");
122 template<
typename TRAITS,
typename LO,
typename GO>
133 std::string fieldName = fieldMap_->find(
scatterFields_[fd].fieldTag().name())->second;
136 subFieldIds_[fd] = rowIndexers_[indexerIds_[fd]]->getFieldNum(fieldName);
144 template<
typename TRAITS,
typename LO,
typename GO>
152 Teuchos::RCP<const BLOC> blockedContainer = Teuchos::rcp_dynamic_cast<
const BLOC>(d.gedc.getDataObject(globalDataKey_));
153 Teuchos::RCP<const ELOC> epetraContainer = Teuchos::rcp_dynamic_cast<
const ELOC>(d.gedc.getDataObject(globalDataKey_));
156 if(blockedContainer!=Teuchos::null)
158 else if(epetraContainer!=Teuchos::null)
159 r_ = Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_f_th());
165 template<
typename TRAITS,
typename LO,
typename GO>
170 using Teuchos::ptrFromRef;
171 using Teuchos::rcp_dynamic_cast;
174 using Thyra::SpmdVectorBase;
178 std::string blockId = this->wda(workset).block_id;
179 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
188 for (std::size_t fieldIndex = 0; fieldIndex <
scatterFields_.size(); fieldIndex++) {
189 int indexerId = indexerIds_[fieldIndex];
190 int subFieldNum = subFieldIds_[fieldIndex];
193 rcp_dynamic_cast<SpmdVectorBase<double> >(r_->getNonconstVectorBlock(indexerId))->getNonconstLocalData(ptrFromRef(local_r));
195 auto subRowIndexer = rowIndexers_[indexerId];
196 const std::vector<int> & elmtOffset = subRowIndexer->getGIDFieldOffsets(blockId,subFieldNum);
199 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
200 std::size_t cellLocalId = localCellIds[worksetCellIndex];
202 const std::vector<LO> & LIDs = subRowIndexer->getElementLIDs(cellLocalId);
218 template<
typename TRAITS,
typename LO,
typename GO>
223 bool useDiscreteAdjoint)
224 : rowIndexers_(rIndexers)
225 , colIndexers_(cIndexers)
226 , globalDataKey_(
"Residual Scatter Container")
228 std::string scatterName = p.
get<std::string>(
"Scatter Name");
233 const std::vector<std::string>& names =
244 for (std::size_t eq = 0; eq < names.size(); ++eq) {
245 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
252 this->addEvaluatedField(*scatterHolder_);
254 if (p.
isType<std::string>(
"Global Data Key"))
255 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
257 this->setName(scatterName+
" Scatter Tangent");
261 template<
typename TRAITS,
typename LO,
typename GO>
272 std::string fieldName = fieldMap_->find(
scatterFields_[fd].fieldTag().name())->second;
275 subFieldIds_[fd] = rowIndexers_[indexerIds_[fd]]->getFieldNum(fieldName);
283 template<
typename TRAITS,
typename LO,
typename GO>
291 Teuchos::RCP<const BLOC> blockedContainer = Teuchos::rcp_dynamic_cast<
const BLOC>(d.gedc.getDataObject(globalDataKey_));
292 Teuchos::RCP<const ELOC> epetraContainer = Teuchos::rcp_dynamic_cast<
const ELOC>(d.gedc.getDataObject(globalDataKey_));
295 if(blockedContainer!=Teuchos::null)
297 else if(epetraContainer!=Teuchos::null)
298 r_ = Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_f_th());
304 template<
typename TRAITS,
typename LO,
typename GO>
311 using Teuchos::ptrFromRef;
312 using Teuchos::rcp_dynamic_cast;
315 using Thyra::SpmdVectorBase;
319 std::string blockId = this->wda(workset).block_id;
320 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
329 for (std::size_t fieldIndex = 0; fieldIndex <
scatterFields_.size(); fieldIndex++) {
330 int indexerId = indexerIds_[fieldIndex];
331 int subFieldNum = subFieldIds_[fieldIndex];
334 rcp_dynamic_cast<SpmdVectorBase<double> >(r_->getNonconstVectorBlock(indexerId))->getNonconstLocalData(ptrFromRef(local_r));
336 auto subRowIndexer = rowIndexers_[indexerId];
337 const std::vector<int> & elmtOffset = subRowIndexer->getGIDFieldOffsets(blockId,subFieldNum);
340 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
341 std::size_t cellLocalId = localCellIds[worksetCellIndex];
343 const std::vector<LO> & LIDs = subRowIndexer->getElementLIDs(cellLocalId);
359 template<
typename TRAITS,
typename LO,
typename GO>
364 bool useDiscreteAdjoint)
365 : rowIndexers_(rIndexers)
366 , colIndexers_(cIndexers)
367 , globalDataKey_(
"Residual Scatter Container")
368 , useDiscreteAdjoint_(useDiscreteAdjoint)
370 std::string scatterName = p.
get<std::string>(
"Scatter Name");
375 const std::vector<std::string>& names =
386 for (std::size_t eq = 0; eq < names.size(); ++eq) {
387 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
394 this->addEvaluatedField(*scatterHolder_);
396 if (p.
isType<std::string>(
"Global Data Key"))
397 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
398 if (p.
isType<
bool>(
"Use Discrete Adjoint"))
399 useDiscreteAdjoint = p.
get<
bool>(
"Use Discrete Adjoint");
402 if(useDiscreteAdjoint)
405 if(colIndexers_.size()==0)
406 colIndexers_ = rowIndexers_;
408 this->setName(scatterName+
" Scatter Residual BlockedEpetra (Jacobian)");
412 template<
typename TRAITS,
typename LO,
typename GO>
423 std::string fieldName = fieldMap_->find(
scatterFields_[fd].fieldTag().name())->second;
426 subFieldIds_[fd] = rowIndexers_[indexerIds_[fd]]->getFieldNum(fieldName);
434 template<
typename TRAITS,
typename LO,
typename GO>
439 using Teuchos::rcp_dynamic_cast;
445 RCP<const BLOC> blockedContainer = rcp_dynamic_cast<
const BLOC>(d.gedc.getDataObject(globalDataKey_));
446 RCP<const ELOC> epetraContainer = rcp_dynamic_cast<
const ELOC>(d.gedc.getDataObject(globalDataKey_));
449 if(blockedContainer!=Teuchos::null) {
453 else if(epetraContainer!=Teuchos::null) {
455 if(epetraContainer->get_f_th()!=Teuchos::null)
456 r_ = Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_f_th());
467 template<
typename TRAITS,
typename LO,
typename GO>
473 using Teuchos::ptrFromRef;
474 using Teuchos::rcp_dynamic_cast;
477 using Thyra::SpmdVectorBase;
481 std::vector<double> jacRow;
484 std::string blockId = this->wda(workset).block_id;
485 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
487 int numFieldBlocks = Teuchos::as<int>(colIndexers_.size());
489 std::vector<int> blockOffsets;
496 for(std::size_t fieldIndex = 0; fieldIndex <
scatterFields_.size(); fieldIndex++) {
497 int rowIndexer = indexerIds_[fieldIndex];
498 int subFieldNum = subFieldIds_[fieldIndex];
501 if(r_!=Teuchos::null)
502 rcp_dynamic_cast<SpmdVectorBase<double> >(r_->getNonconstVectorBlock(rowIndexer))->getNonconstLocalData(ptrFromRef(local_r));
504 auto subRowIndexer = rowIndexers_[rowIndexer];
505 const std::vector<int> & elmtOffset = subRowIndexer->getGIDFieldOffsets(blockId,subFieldNum);
508 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
509 std::size_t cellLocalId = localCellIds[worksetCellIndex];
511 const std::vector<LO> & rLIDs = subRowIndexer->getElementLIDs(cellLocalId);
514 for(std::size_t rowBasisNum = 0; rowBasisNum < elmtOffset.size(); rowBasisNum++) {
516 int rowOffset = elmtOffset[rowBasisNum];
517 int r_lid = rLIDs[rowOffset];
520 if(local_r!=Teuchos::null)
521 local_r[r_lid] += (scatterField.val());
524 jacRow.
resize(scatterField.size());
527 if(scatterField.size() == 0)
530 for(
int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
531 jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex);
534 for(
int colIndexer=0;colIndexer<numFieldBlocks;colIndexer++) {
535 int start = blockOffsets[colIndexer];
536 int end = blockOffsets[colIndexer+1];
541 auto subColIndexer = colIndexers_[colIndexer];
542 const std::vector<LO> & cLIDs = subColIndexer->getElementLIDs(cellLocalId);
547 std::pair<int,int> blockIndex = useDiscreteAdjoint_ ? std::make_pair(colIndexer,rowIndexer)
548 : std::make_pair(rowIndexer,colIndexer);
552 if(subJac==Teuchos::null) {
560 subJac = rcp_dynamic_cast<Epetra_CrsMatrix>(eOp,
true);
561 jacEpetraBlocks[blockIndex] = subJac;
565 if(!useDiscreteAdjoint_) {
566 int err = subJac->SumIntoMyValues(r_lid, end-start, &jacRow[start],&cLIDs[0]);
569 std::stringstream ss;
570 ss <<
"Failed inserting row: " <<
"LID = " << r_lid <<
"): ";
571 for(
int i=start;i<end;i++)
572 ss << cLIDs[i] <<
" ";
574 ss <<
"Into block " << rowIndexer <<
", " << colIndexer << std::endl;
576 ss <<
"scatter field = ";
584 for(std::size_t c=0;c<cLIDs.size();c++) {
585 int err = subJac->SumIntoMyValues(cLIDs[c], 1, &jacRow[start+c],&r_lid);
ScatterResidual_BlockedEpetra(const Teuchos::ParameterList &p)
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_
T & get(const std::string &name, T def_value)
bool is_null(const std::shared_ptr< T > &p)
void computeBlockOffsets(const std::string &blockId, const std::vector< Teuchos::RCP< UniqueGlobalIndexer< LocalOrdinalT, GlobalOrdinalT > > > &ugis, std::vector< int > &blockOffsets)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
PHX::MDField< ScalarT > vector
void evaluateFields(typename TRAITS::EvalData d)
void resize(const size_type n, const T &val=T())
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
bool isType(const std::string &name) const
Pushes residual values into the residual vector for a Newton-based solve.
panzer::Traits::Jacobian::ScalarT ScalarT
Teuchos::RCP< const panzer::PureBasis > basis
Interpolates basis DOF values to IP DOF values.
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &vm)
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_ASSERT_EQUALITY(val1, val2)