43#ifndef PANZER_SCATTER_RESIDUAL_BLOCKEDEPETRA_IMPL_HPP
44#define PANZER_SCATTER_RESIDUAL_BLOCKEDEPETRA_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"
63#include "Thyra_SpmdVectorBase.hpp"
64#include "Thyra_ProductVectorBase.hpp"
65#include "Thyra_DefaultProductVector.hpp"
66#include "Thyra_BlockedLinearOpBase.hpp"
67#include "Thyra_DefaultBlockedLinearOp.hpp"
68#include "Thyra_get_Epetra_Operator.hpp"
70#include "Phalanx_DataLayout_MDALayout.hpp"
72#include "Teuchos_FancyOStream.hpp"
78template<
typename TRAITS,
typename LO,
typename GO>
81 const std::vector<Teuchos::RCP<const GlobalIndexer> > & cIndexers,
82 const Teuchos::ParameterList& p,
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");
99 Teuchos::RCP<PHX::DataLayout> dl =
100 p.get< Teuchos::RCP<const panzer::PureBasis> >(
"Basis")->functional;
103 scatterFields_.resize(names.size());
104 for (std::size_t eq = 0; eq < names.size(); ++eq) {
105 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
108 this->addDependentField(scatterFields_[eq]);
112 this->addEvaluatedField(*scatterHolder_);
114 if (p.isType<std::string>(
"Global Data Key"))
115 globalDataKey_ = p.get<std::string>(
"Global Data Key");
117 this->setName(scatterName+
" Scatter Residual");
121template<
typename TRAITS,
typename LO,
typename GO>
126 indexerIds_.resize(scatterFields_.size());
127 subFieldIds_.resize(scatterFields_.size());
130 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
132 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
135 subFieldIds_[fd] = rowIndexers_[indexerIds_[fd]]->getFieldNum(fieldName);
140template<
typename TRAITS,
typename LO,
typename GO>
148 Teuchos::RCP<const BLOC> blockedContainer = Teuchos::rcp_dynamic_cast<const BLOC>(d.gedc->getDataObject(globalDataKey_));
149 Teuchos::RCP<const ELOC> epetraContainer = Teuchos::rcp_dynamic_cast<const ELOC>(d.gedc->getDataObject(globalDataKey_));
152 if(blockedContainer!=Teuchos::null)
153 r_ = Teuchos::rcp_dynamic_cast<Thyra::ProductVectorBase<double> >(blockedContainer->get_f());
154 else if(epetraContainer!=Teuchos::null)
155 r_ = Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_f_th());
157 TEUCHOS_ASSERT(r_!=Teuchos::null);
161template<
typename TRAITS,
typename LO,
typename GO>
166 using Teuchos::ptrFromRef;
167 using Teuchos::rcp_dynamic_cast;
170 using Thyra::SpmdVectorBase;
174 std::string blockId = this->wda(workset).block_id;
175 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
183 Teuchos::ArrayRCP<double> local_r;
184 for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
185 int indexerId = indexerIds_[fieldIndex];
186 int subFieldNum = subFieldIds_[fieldIndex];
189 rcp_dynamic_cast<SpmdVectorBase<double> >(r_->getNonconstVectorBlock(indexerId))->getNonconstLocalData(ptrFromRef(local_r));
191 auto subRowIndexer = rowIndexers_[indexerId];
192 const std::vector<int> & elmtOffset = subRowIndexer->getGIDFieldOffsets(blockId,subFieldNum);
194 auto field = PHX::as_view(scatterFields_[fieldIndex]);
195 auto field_h = Kokkos::create_mirror_view(
field);
196 Kokkos::deep_copy(field_h,
field);
199 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
200 std::size_t cellLocalId = localCellIds[worksetCellIndex];
202 auto LIDs = subRowIndexer->getElementLIDs(cellLocalId);
203 auto LIDs_h = Kokkos::create_mirror_view(LIDs);
204 Kokkos::deep_copy(LIDs_h, LIDs);
207 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
208 int offset = elmtOffset[basis];
209 int lid = LIDs_h[offset];
210 local_r[lid] += field_h(worksetCellIndex,basis);
220template<
typename TRAITS,
typename LO,
typename GO>
223 const std::vector<Teuchos::RCP<const GlobalIndexer> > & cIndexers,
224 const Teuchos::ParameterList& p,
226 : rowIndexers_(rIndexers)
227 , colIndexers_(cIndexers)
228 , globalDataKey_(
"Residual Scatter Container")
230 std::string scatterName = p.get<std::string>(
"Scatter Name");
232 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
235 const std::vector<std::string>& names =
236 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
239 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
241 Teuchos::RCP<PHX::DataLayout> dl =
242 p.get< Teuchos::RCP<const panzer::PureBasis> >(
"Basis")->functional;
245 scatterFields_.resize(names.size());
246 for (std::size_t eq = 0; eq < names.size(); ++eq) {
247 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
250 this->addDependentField(scatterFields_[eq]);
254 this->addEvaluatedField(*scatterHolder_);
256 if (p.isType<std::string>(
"Global Data Key"))
257 globalDataKey_ = p.get<std::string>(
"Global Data Key");
259 this->setName(scatterName+
" Scatter Tangent");
263template<
typename TRAITS,
typename LO,
typename GO>
268 indexerIds_.resize(scatterFields_.size());
269 subFieldIds_.resize(scatterFields_.size());
272 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
274 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
277 subFieldIds_[fd] = rowIndexers_[indexerIds_[fd]]->getFieldNum(fieldName);
282template<
typename TRAITS,
typename LO,
typename GO>
290 Teuchos::RCP<const BLOC> blockedContainer = Teuchos::rcp_dynamic_cast<const BLOC>(d.gedc->getDataObject(globalDataKey_));
291 Teuchos::RCP<const ELOC> epetraContainer = Teuchos::rcp_dynamic_cast<const ELOC>(d.gedc->getDataObject(globalDataKey_));
294 if(blockedContainer!=Teuchos::null)
295 r_ = Teuchos::rcp_dynamic_cast<Thyra::ProductVectorBase<double> >(blockedContainer->get_f());
296 else if(epetraContainer!=Teuchos::null)
297 r_ = Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_f_th());
299 TEUCHOS_ASSERT(r_!=Teuchos::null);
303template<
typename TRAITS,
typename LO,
typename GO>
307 TEUCHOS_ASSERT(
false);
310 using Teuchos::ptrFromRef;
311 using Teuchos::rcp_dynamic_cast;
314 using Thyra::SpmdVectorBase;
318 std::string blockId = this->wda(workset).block_id;
319 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
327 Teuchos::ArrayRCP<double> local_r;
328 for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
329 int indexerId = indexerIds_[fieldIndex];
330 int subFieldNum = subFieldIds_[fieldIndex];
333 rcp_dynamic_cast<SpmdVectorBase<double> >(r_->getNonconstVectorBlock(indexerId))->getNonconstLocalData(ptrFromRef(local_r));
335 auto subRowIndexer = rowIndexers_[indexerId];
336 const std::vector<int> & elmtOffset = subRowIndexer->getGIDFieldOffsets(blockId,subFieldNum);
339 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
340 std::size_t cellLocalId = localCellIds[worksetCellIndex];
342 auto LIDs = subRowIndexer->getElementLIDs(cellLocalId);
345 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
346 int offset = elmtOffset[basis];
347 int lid = LIDs[offset];
348 local_r[lid] += (scatterFields_[fieldIndex])(worksetCellIndex,basis).val();
358template<
typename TRAITS,
typename LO,
typename GO>
361 const std::vector<Teuchos::RCP<const GlobalIndexer> > & cIndexers,
362 const Teuchos::ParameterList& p,
363 bool useDiscreteAdjoint)
364 : rowIndexers_(rIndexers)
365 , colIndexers_(cIndexers)
366 , globalDataKey_(
"Residual Scatter Container")
367 , useDiscreteAdjoint_(useDiscreteAdjoint)
369 std::string scatterName = p.get<std::string>(
"Scatter Name");
371 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
374 const std::vector<std::string>& names =
375 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
378 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
380 Teuchos::RCP<PHX::DataLayout> dl =
381 p.get< Teuchos::RCP<const panzer::PureBasis> >(
"Basis")->functional;
384 scatterFields_.resize(names.size());
385 for (std::size_t eq = 0; eq < names.size(); ++eq) {
386 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
389 this->addDependentField(scatterFields_[eq]);
393 this->addEvaluatedField(*scatterHolder_);
395 if (p.isType<std::string>(
"Global Data Key"))
396 globalDataKey_ = p.get<std::string>(
"Global Data Key");
397 if (p.isType<
bool>(
"Use Discrete Adjoint"))
398 useDiscreteAdjoint = p.get<
bool>(
"Use Discrete Adjoint");
401 if(useDiscreteAdjoint)
402 { TEUCHOS_ASSERT(colIndexers_.size()==0); }
404 if(colIndexers_.size()==0)
405 colIndexers_ = rowIndexers_;
407 this->setName(scatterName+
" Scatter Residual BlockedEpetra (Jacobian)");
411template<
typename TRAITS,
typename LO,
typename GO>
416 indexerIds_.resize(scatterFields_.size());
417 subFieldIds_.resize(scatterFields_.size());
420 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
422 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
425 subFieldIds_[fd] = rowIndexers_[indexerIds_[fd]]->getFieldNum(fieldName);
430template<
typename TRAITS,
typename LO,
typename GO>
435 using Teuchos::rcp_dynamic_cast;
441 RCP<const BLOC> blockedContainer = rcp_dynamic_cast<const BLOC>(d.gedc->getDataObject(globalDataKey_));
442 RCP<const ELOC> epetraContainer = rcp_dynamic_cast<const ELOC>(d.gedc->getDataObject(globalDataKey_));
445 if(blockedContainer!=Teuchos::null) {
446 r_ = Teuchos::rcp_dynamic_cast<Thyra::ProductVectorBase<double> >(blockedContainer->get_f());
447 Jac_ = rcp_dynamic_cast<Thyra::BlockedLinearOpBase<double> >(blockedContainer->get_A());
449 else if(epetraContainer!=Teuchos::null) {
451 if(epetraContainer->get_f_th()!=Teuchos::null)
452 r_ = Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_f_th());
455 RCP<Thyra::LinearOpBase<double> > J = blockedContainer->get_A_th();
456 Jac_ = rcp_dynamic_cast<Thyra::BlockedLinearOpBase<double> >(Thyra::nonconstBlock1x1(J));
459 TEUCHOS_ASSERT(Jac_!=Teuchos::null);
463template<
typename TRAITS,
typename LO,
typename GO>
468 using Teuchos::ArrayRCP;
469 using Teuchos::ptrFromRef;
470 using Teuchos::rcp_dynamic_cast;
473 using Thyra::SpmdVectorBase;
477 std::vector<double> jacRow;
480 std::string blockId = this->wda(workset).block_id;
481 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
483 int numFieldBlocks = Teuchos::as<int>(colIndexers_.size());
485 std::vector<int> blockOffsets;
488 std::unordered_map<std::pair<int,int>,Teuchos::RCP<Epetra_CrsMatrix>,
panzer::pair_hash> jacEpetraBlocks;
491 Teuchos::ArrayRCP<double> local_r;
492 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
493 int rowIndexer = indexerIds_[fieldIndex];
494 int subFieldNum = subFieldIds_[fieldIndex];
497 if(r_!=Teuchos::null)
498 rcp_dynamic_cast<SpmdVectorBase<double> >(r_->getNonconstVectorBlock(rowIndexer))->getNonconstLocalData(ptrFromRef(local_r));
500 auto subRowIndexer = rowIndexers_[rowIndexer];
501 const std::vector<int> & elmtOffset = subRowIndexer->getGIDFieldOffsets(blockId,subFieldNum);
503 auto field = scatterFields_[fieldIndex].get_view();
504 auto field_h = Kokkos::create_mirror_view(
field);
505 Kokkos::deep_copy(field_h,
field);
507 auto rLIDs = subRowIndexer->getLIDs();
508 auto rLIDs_h = Kokkos::create_mirror_view(rLIDs);
509 Kokkos::deep_copy(rLIDs_h, rLIDs);
512 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
513 std::size_t cellLocalId = localCellIds[worksetCellIndex];
516 for(std::size_t rowBasisNum = 0; rowBasisNum < elmtOffset.size(); rowBasisNum++) {
517 const ScalarT scatterField = field_h(worksetCellIndex,rowBasisNum);
518 int rowOffset = elmtOffset[rowBasisNum];
519 int r_lid = rLIDs_h(cellLocalId, rowOffset);
522 if(local_r!=Teuchos::null)
523 local_r[r_lid] += (scatterField.val());
526 jacRow.resize(scatterField.size());
529 if(scatterField.size() == 0)
532 for(
int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
533 jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex);
536 for(
int colIndexer=0;colIndexer<numFieldBlocks;colIndexer++) {
537 int start = blockOffsets[colIndexer];
538 int end = blockOffsets[colIndexer+1];
543 auto subColIndexer = colIndexers_[colIndexer];
544 auto cLIDs = subColIndexer->getElementLIDs(cellLocalId);
545 auto cLIDs_h = Kokkos::create_mirror_view(cLIDs);
546 Kokkos::deep_copy(cLIDs_h, cLIDs);
548 TEUCHOS_ASSERT(end-start==Teuchos::as<int>(cLIDs.size()));
551 std::pair<int,int> blockIndex = useDiscreteAdjoint_ ? std::make_pair(colIndexer,rowIndexer)
552 : std::make_pair(rowIndexer,colIndexer);
553 Teuchos::RCP<Epetra_CrsMatrix> subJac = jacEpetraBlocks[blockIndex];
556 if(subJac==Teuchos::null) {
557 Teuchos::RCP<Thyra::LinearOpBase<double> > tOp = Jac_->getNonconstBlock(blockIndex.first,blockIndex.second);
560 if(Teuchos::is_null(tOp))
563 Teuchos::RCP<Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*tOp);
564 subJac = rcp_dynamic_cast<Epetra_CrsMatrix>(eOp,
true);
565 jacEpetraBlocks[blockIndex] = subJac;
569 if(!useDiscreteAdjoint_) {
570 int err = subJac->SumIntoMyValues(r_lid, end-start, &jacRow[start],&cLIDs_h[0]);
573 std::stringstream ss;
574 ss <<
"Failed inserting row: " <<
"LID = " << r_lid <<
"): ";
575 for(
int i=start;i<end;i++)
576 ss << cLIDs_h[i] <<
" ";
578 ss <<
"Into block " << rowIndexer <<
", " << colIndexer << std::endl;
580 ss <<
"scatter field = ";
581 scatterFields_[fieldIndex].print(ss);
584 TEUCHOS_TEST_FOR_EXCEPTION(err!=0,std::runtime_error,ss.str());
588 for(std::size_t c=0;c<cLIDs.size();c++) {
589 int err = subJac->SumIntoMyValues(cLIDs_h[c], 1, &jacRow[start+c],&r_lid);
590 TEUCHOS_ASSERT_EQUALITY(err,0);
PHX::MDField< ScalarT, panzer::Cell, panzer::BASIS > field
A field to which we'll contribute, or in which we'll store, the result of computing this integral.
panzer::Traits::Jacobian::ScalarT ScalarT
Pushes residual values into the residual vector for a Newton-based solve.
ScatterResidual_BlockedEpetra(const Teuchos::ParameterList &p)
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &vm)
void evaluateFields(typename TRAITS::EvalData d)
void computeBlockOffsets(const std::string &blockId, const std::vector< Teuchos::RCP< GlobalIndexer > > &ugis, std::vector< int > &blockOffsets)
int getFieldBlock(const std::string &fieldName, const std::vector< Teuchos::RCP< const GlobalIndexer > > &ugis)