43 #ifndef __Panzer_GatherSolution_BlockedEpetra_impl_hpp__
44 #define __Panzer_GatherSolution_BlockedEpetra_impl_hpp__
62 #include "Phalanx_DataLayout.hpp"
65 #include "Teuchos_Assert.hpp"
66 #include "Teuchos_FancyOStream.hpp"
69 #include "Thyra_ProductVectorBase.hpp"
70 #include "Thyra_SpmdVectorBase.hpp"
77 template<
typename TRAITS,
typename LO,
typename GO>
86 hasTangentFields_(false)
95 using vvstring = std::vector<std::vector<std::string>>;
111 MDField<ScalarT, Cell, NODE>(names[fd], basis->functional);
112 this->addEvaluatedField(gatherFields_[fd]);
116 if (tangentFieldNames.size() > 0)
119 hasTangentFields_ =
true;
120 tangentFields_.resize(numFields);
123 int numTangentFields(tangentFieldNames[fd].size());
124 tangentFields_[fd].resize(numTangentFields);
125 for (
int i(0); i < numTangentFields; ++i)
127 tangentFields_[fd][i] =
128 MDField<const ScalarT, Cell, NODE>(tangentFieldNames[fd][i],
130 this->addDependentField(tangentFields_[fd][i]);
136 string firstName(
"<none>");
138 firstName = names[0];
139 string n(
"GatherSolution (BlockedEpetra): " + firstName +
" (" +
140 print<EvalT>() +
")");
149 template<
typename TRAITS,
typename LO,
typename GO>
154 typename TRAITS::SetupData ,
166 const string& fieldName(indexerNames_[fd]);
168 subFieldIds_[fd] = indexers_[indexerIds_[fd]]->getFieldNum(fieldName);
171 indexerNames_.clear();
179 template<
typename TRAITS,
typename LO,
typename GO>
184 typename TRAITS::PreEvalData d)
186 using std::logic_error;
189 using Teuchos::rcp_dynamic_cast;
198 string post(useTimeDerivativeSolutionVector_ ?
" - Xdot" :
" - X");
199 if (d.gedc->containsDataObject(globalDataKey_ + post))
201 ged = d.gedc->getDataObject(globalDataKey_ + post);
202 xBvRoGed_ = rcp_dynamic_cast<BVROGED>(ged,
true);
208 ged = d.gedc->getDataObject(globalDataKey_);
211 auto roGed = rcp_dynamic_cast<
const BVROGED>(ged);
212 auto beLoc = rcp_dynamic_cast<
const BELOC>(ged);
213 if (not roGed.is_null())
214 xBvRoGed_ = rcp_dynamic_cast<BVROGED>(ged,
true);
215 else if (not beLoc.is_null())
217 if (useTimeDerivativeSolutionVector_)
218 x_ = rcp_dynamic_cast<ProductVectorBase<double>>(beLoc->get_dxdt());
220 x_ = rcp_dynamic_cast<ProductVectorBase<double>>(beLoc->get_x());
222 "Residual: Can't find the x_ vector in GEDC \"" << globalDataKey_ <<
223 "\" (" << post <<
"). A cast failed for " << ged <<
". Type is " <<
234 template<
typename TRAITS,
typename LO,
typename GO>
239 typename TRAITS::EvalData workset)
246 using Teuchos::ptrFromRef;
248 using Teuchos::rcp_dynamic_cast;
250 using Thyra::SpmdVectorBase;
254 string blockId(this->wda(workset).block_id);
255 const vector<size_t>& localCellIds = this->wda(workset).cell_local_ids;
256 int numFields(gatherFields_.size()), numCells(localCellIds.size());
259 for (
int fieldInd(0); fieldInd <
numFields; ++fieldInd)
261 MDField<ScalarT, Cell, NODE>&
field = gatherFields_[fieldInd];
262 auto field_h = Kokkos::create_mirror_view(field.get_static_view());
264 int indexerId(indexerIds_[fieldInd]),
265 subFieldNum(subFieldIds_[fieldInd]);
271 if(not x_.is_null()) {
272 rcp_dynamic_cast<SpmdVectorBase<double>>(x_->
273 getNonconstVectorBlock(indexerId))->getLocalData(ptrFromRef(x));
276 xEvRoGed = xBvRoGed_->getGEDBlock(indexerId);
279 auto subRowIndexer = indexers_[indexerId];
280 const vector<int>& elmtOffset =
281 subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
282 int numBases(elmtOffset.size());
284 auto LIDs = subRowIndexer->getLIDs();
285 auto LIDs_h = Kokkos::create_mirror_view(LIDs);
286 Kokkos::deep_copy(LIDs_h, LIDs);
289 for (
int cell(0); cell < numCells; ++cell)
291 LO cellLocalId = localCellIds[cell];
294 for (
int basis(0); basis < numBases; ++basis)
296 int offset(elmtOffset[basis]), lid(LIDs_h(cellLocalId, offset));
298 field_h(cell, basis) = (*xEvRoGed)[lid];
300 field_h(cell, basis) = x[lid];
303 Kokkos::deep_copy(field.get_static_view(), field_h);
312 template<
typename TRAITS,
typename LO,
typename GO>
320 hasTangentFields_(false)
329 using vvstring = std::vector<std::vector<std::string>>;
345 MDField<ScalarT, Cell, NODE>(names[fd], basis->functional);
346 this->addEvaluatedField(gatherFields_[fd]);
350 if (tangentFieldNames.size() > 0)
353 hasTangentFields_ =
true;
354 tangentFields_.resize(numFields);
357 int numTangentFields(tangentFieldNames[fd].size());
358 tangentFields_[fd].resize(numTangentFields);
359 for (
int i(0); i < numTangentFields; ++i)
361 tangentFields_[fd][i] =
362 MDField<const ScalarT, Cell, NODE>(tangentFieldNames[fd][i],
364 this->addDependentField(tangentFields_[fd][i]);
370 string firstName(
"<none>");
372 firstName = names[0];
373 string n(
"GatherSolution Tangent (BlockedEpetra): " + firstName +
" (" +
374 print<EvalT>() +
")");
383 template<
typename TRAITS,
typename LO,
typename GO>
387 typename TRAITS::SetupData ,
399 const string& fieldName(indexerNames_[fd]);
401 subFieldIds_[fd] = indexers_[indexerIds_[fd]]->getFieldNum(fieldName);
404 indexerNames_.clear();
412 template<
typename TRAITS,
typename LO,
typename GO>
416 typename TRAITS::PreEvalData d)
418 using std::logic_error;
421 using Teuchos::rcp_dynamic_cast;
430 string post(useTimeDerivativeSolutionVector_ ?
" - Xdot" :
" - X");
431 if (d.gedc->containsDataObject(globalDataKey_ + post))
433 ged = d.gedc->getDataObject(globalDataKey_ + post);
434 xBvRoGed_ = rcp_dynamic_cast<BVROGED>(ged,
true);
440 ged = d.gedc->getDataObject(globalDataKey_);
443 auto roGed = rcp_dynamic_cast<
const BVROGED>(ged);
444 auto beLoc = rcp_dynamic_cast<
const BELOC>(ged);
445 if (not roGed.is_null())
446 xBvRoGed_ = rcp_dynamic_cast<BVROGED>(ged,
true);
447 else if (not beLoc.is_null())
449 if (useTimeDerivativeSolutionVector_)
450 x_ = rcp_dynamic_cast<ProductVectorBase<double>>(beLoc->get_dxdt());
452 x_ = rcp_dynamic_cast<ProductVectorBase<double>>(beLoc->get_x());
454 "Tangent: Can't find the x_ vector in GEDC \"" << globalDataKey_ <<
455 "\" (" << post <<
"). A cast failed for " << ged <<
". Type is " <<
466 template<
typename TRAITS,
typename LO,
typename GO>
470 typename TRAITS::EvalData workset)
478 using Teuchos::ptrFromRef;
480 using Teuchos::rcp_dynamic_cast;
482 using Thyra::SpmdVectorBase;
486 vector<pair<int, GO>> GIDs;
487 string blockId(this->wda(workset).block_id);
488 const vector<size_t>& localCellIds = this->wda(workset).cell_local_ids;
489 int numFields(gatherFields_.size()), numCells(localCellIds.size());
494 for (
int fieldInd(0); fieldInd <
numFields; ++fieldInd)
496 MDField<ScalarT, Cell, NODE>&
field = gatherFields_[fieldInd];
497 int indexerId(indexerIds_[fieldInd]),
498 subFieldNum(subFieldIds_[fieldInd]);
501 auto xEvRoGed = xBvRoGed_->getGEDBlock(indexerId);
502 auto subRowIndexer = indexers_[indexerId];
503 const vector<int>& elmtOffset =
504 subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
505 int numBases(elmtOffset.size());
508 for (
int cell(0); cell < numCells; ++cell)
510 LO cellLocalId = localCellIds[cell];
511 auto LIDs = subRowIndexer->getElementLIDs(cellLocalId);
514 for (
int basis(0); basis < numBases; ++basis)
516 int offset(elmtOffset[basis]), lid(LIDs[offset]);
517 field(cell, basis) = (*xEvRoGed)[lid];
525 for (
int fieldInd(0); fieldInd <
numFields; ++fieldInd)
527 MDField<ScalarT, Cell, NODE>&
field = gatherFields_[fieldInd];
528 int indexerId(indexerIds_[fieldInd]),
529 subFieldNum(subFieldIds_[fieldInd]);
533 rcp_dynamic_cast<SpmdVectorBase<double>>(x_->
534 getNonconstVectorBlock(indexerId))->getLocalData(ptrFromRef(x));
535 auto subRowIndexer = indexers_[indexerId];
536 const vector<int>& elmtOffset =
537 subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
538 int numBases(elmtOffset.size());
541 for (
int cell(0); cell < numCells; ++cell)
543 LO cellLocalId = localCellIds[cell];
544 auto LIDs = subRowIndexer->getElementLIDs(cellLocalId);
547 for (
int basis(0); basis < numBases; ++basis)
549 int offset(elmtOffset[basis]), lid(LIDs[offset]);
550 field(cell, basis) = x[lid];
557 if (hasTangentFields_)
560 for (
int fieldInd(0); fieldInd <
numFields; ++fieldInd)
562 MDField<ScalarT, Cell, NODE>&
field = gatherFields_[fieldInd];
563 int indexerId(indexerIds_[fieldInd]),
564 subFieldNum(subFieldIds_[fieldInd]);
565 auto subRowIndexer = indexers_[indexerId];
566 const vector<int>& elmtOffset =
567 subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
568 int numBases(elmtOffset.size());
571 for (
int cell(0); cell < numCells; ++cell)
574 for (
int basis(0); basis < numBases; ++basis)
576 int numTangentFields(tangentFields_[fieldInd].size());
577 for (
int i(0); i < numTangentFields; ++i)
578 field(cell, basis).fastAccessDx(i) =
579 tangentFields_[fieldInd][i](cell, basis).val();
591 template<
typename TRAITS,
typename LO,
typename GO>
624 MDField<ScalarT, Cell, NODE> f(names[fd], basis->functional);
625 gatherFields_[fd] = f;
626 this->addEvaluatedField(gatherFields_[fd]);
630 string firstName(
"<none>"), n(
"GatherSolution (BlockedEpetra");
632 firstName = names[0];
633 if (disableSensitivities_)
634 n +=
", No Sensitivities";
635 n +=
"): " + firstName +
" (" + print<EvalT>() +
")";
644 template<
typename TRAITS,
typename LO,
typename GO>
649 typename TRAITS::SetupData ,
661 const string& fieldName(indexerNames_[fd]);
663 subFieldIds_[fd] = indexers_[indexerIds_[fd]]->getFieldNum(fieldName);
666 indexerNames_.clear();
674 template<
typename TRAITS,
typename LO,
typename GO>
679 typename TRAITS::PreEvalData d)
682 using std::logic_error;
685 using Teuchos::rcp_dynamic_cast;
693 applySensitivities_ =
false;
694 if ((not disableSensitivities_ ) and
695 (d.first_sensitivities_name == sensitivitiesName_))
696 applySensitivities_ =
true;
700 string post(useTimeDerivativeSolutionVector_ ?
" - Xdot" :
" - X");
701 if (d.gedc->containsDataObject(globalDataKey_ + post))
703 ged = d.gedc->getDataObject(globalDataKey_ + post);
704 xBvRoGed_ = rcp_dynamic_cast<BVROGED>(ged,
true);
710 ged = d.gedc->getDataObject(globalDataKey_);
713 auto roGed = rcp_dynamic_cast<
const BVROGED>(ged);
714 auto beLoc = rcp_dynamic_cast<
const BELOC>(ged);
715 if (not roGed.is_null())
716 xBvRoGed_ = rcp_dynamic_cast<BVROGED>(ged,
true);
717 else if (not beLoc.is_null())
719 if (useTimeDerivativeSolutionVector_)
720 x_ = rcp_dynamic_cast<ProductVectorBase<double>>(beLoc->get_dxdt());
722 x_ = rcp_dynamic_cast<ProductVectorBase<double>>(beLoc->get_x());
724 "Jacobian: Can't find x vector in GEDC \"" << globalDataKey_ <<
725 "\" (" << post <<
"). A cast failed for " << ged <<
". Type is " <<
736 template<
typename TRAITS,
typename LO,
typename GO>
741 typename TRAITS::EvalData workset)
748 using Teuchos::ptrFromRef;
750 using Teuchos::rcp_dynamic_cast;
752 using Thyra::SpmdVectorBase;
756 string blockId(this->wda(workset).block_id);
757 const vector<size_t>& localCellIds = this->wda(workset).cell_local_ids;
758 int numFields(gatherFields_.size()), numCells(localCellIds.size());
760 if (applySensitivities_)
762 if ((useTimeDerivativeSolutionVector_) and (gatherSeedIndex_ < 0))
763 seedValue = workset.alpha;
764 else if (gatherSeedIndex_ < 0)
765 seedValue = workset.beta;
766 else if (not useTimeDerivativeSolutionVector_)
767 seedValue = workset.gather_seeds[gatherSeedIndex_];
775 if (not applySensitivities_)
778 vector<int> blockOffsets;
780 int numDerivs(blockOffsets[blockOffsets.size() - 1]);
787 for (
int fieldInd(0); fieldInd <
numFields; ++fieldInd)
789 MDField<ScalarT, Cell, NODE>&
field = gatherFields_[fieldInd];
790 auto field_h = Kokkos::create_mirror_view(field.get_view());
792 int indexerId(indexerIds_[fieldInd]), subFieldNum(subFieldIds_[fieldInd]);
797 if(not x_.is_null()) {
798 rcp_dynamic_cast<SpmdVectorBase<double>>(x_->getNonconstVectorBlock(indexerId))->getLocalData(ptrFromRef(x));
800 xEvRoGed = xBvRoGed_->getGEDBlock(indexerId);
803 auto subRowIndexer = indexers_[indexerId];
804 const vector<int>& elmtOffset =
805 subRowIndexer->getGIDFieldOffsets(blockId, subFieldNum);
806 int startBlkOffset(blockOffsets[indexerId]), numBases(elmtOffset.size());
808 auto LIDs = subRowIndexer->getLIDs();
809 auto LIDs_h = Kokkos::create_mirror_view(LIDs);
810 Kokkos::deep_copy(LIDs_h, LIDs);
813 for (
int cell(0); cell < numCells; ++cell)
815 LO cellLocalId = localCellIds[cell];
818 for (
int basis(0); basis < numBases; ++basis)
821 int offset(elmtOffset[basis]), lid(LIDs_h(cellLocalId, offset));
823 field_h(cell, basis) =
ScalarT(numDerivs, (*xEvRoGed)[lid]);
825 field_h(cell, basis) =
ScalarT(numDerivs, x[lid]);
827 field_h(cell, basis).fastAccessDx(startBlkOffset + offset) = seedValue;
830 Kokkos::deep_copy(field.get_static_view(), field_h);
834 #endif // __Panzer_GatherSolution_BlockedEpetra_impl_hpp__
#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)
This class encapsulates the needs of a gather operation to do a halo exchange for blocked vectors...
int getFieldBlock(const std::string &fieldName, const std::vector< Teuchos::RCP< const GlobalIndexer >> &ugis)
void postRegistrationSetup(typename TRAITS::SetupData d, PHX::FieldManager< TRAITS > &fm)
Post-Registration Setup.
GatherSolution_BlockedEpetra(const Teuchos::ParameterList &p)
Constructor.
virtual void preEvaluate(typename Traits::PreEvalData d)=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...
void evaluateFields(typename TRAITS::EvalData d)
Evaluate Fields.
Description and data layouts associated with a particular basis.
#define TEUCHOS_ASSERT(assertion_test)
panzer::Traits::Jacobian::ScalarT ScalarT
The scalar type.
std::string typeName(const T &t)