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"
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"
78 template<
typename TRAITS,
typename LO,
typename GO>
84 : rowIndexers_(rIndexers)
85 , colIndexers_(cIndexers)
86 , globalDataKey_(
"Residual Scatter Container")
88 std::string scatterName = p.
get<std::string>(
"Scatter Name");
93 const std::vector<std::string>& names =
103 scatterFields_.resize(names.size());
104 for (std::size_t eq = 0; eq < names.size(); ++eq) {
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");
121 template<
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);
140 template<
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)
154 else if(epetraContainer!=Teuchos::null)
155 r_ = Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_f_th());
161 template<
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;
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);
220 template<
typename TRAITS,
typename LO,
typename GO>
226 : rowIndexers_(rIndexers)
227 , colIndexers_(cIndexers)
228 , globalDataKey_(
"Residual Scatter Container")
230 std::string scatterName = p.
get<std::string>(
"Scatter Name");
235 const std::vector<std::string>& names =
245 scatterFields_.resize(names.size());
246 for (std::size_t eq = 0; eq < names.size(); ++eq) {
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");
263 template<
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);
282 template<
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)
296 else if(epetraContainer!=Teuchos::null)
297 r_ = Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_f_th());
303 template<
typename TRAITS,
typename LO,
typename GO>
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;
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();
358 template<
typename TRAITS,
typename LO,
typename GO>
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");
374 const std::vector<std::string>& names =
384 scatterFields_.resize(names.size());
385 for (std::size_t eq = 0; eq < names.size(); ++eq) {
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)
404 if(colIndexers_.size()==0)
405 colIndexers_ = rowIndexers_;
407 this->setName(scatterName+
" Scatter Residual BlockedEpetra (Jacobian)");
411 template<
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);
430 template<
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) {
449 else if(epetraContainer!=Teuchos::null) {
451 if(epetraContainer->get_f_th()!=Teuchos::null)
452 r_ = Thyra::castOrCreateNonconstProductVectorBase<double>(epetraContainer->get_f_th());
463 template<
typename TRAITS,
typename LO,
typename GO>
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;
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);
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) {
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);
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);
ScatterResidual_BlockedEpetra(const Teuchos::ParameterList &p)
T & get(const std::string &name, T def_value)
bool is_null(const std::shared_ptr< T > &p)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void computeBlockOffsets(const std::string &blockId, const std::vector< Teuchos::RCP< GlobalIndexer >> &ugis, std::vector< int > &blockOffsets)
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)
int getFieldBlock(const std::string &fieldName, const std::vector< Teuchos::RCP< const GlobalIndexer >> &ugis)
int SumIntoMyValues(int MyRow, int NumEntries, const double *Values, const int *Indices)
virtual void preEvaluate(typename Traits::PreEvalData d)=0
panzer::Traits::Jacobian::ScalarT ScalarT
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...
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &vm)
bool isType(const std::string &name) const
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_ASSERT_EQUALITY(val1, val2)