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" 59 #include "Thyra_ProductVectorBase.hpp" 60 #include "Thyra_BlockedLinearOpBase.hpp" 61 #include "Thyra_TpetraVector.hpp" 62 #include "Thyra_TpetraLinearOp.hpp" 63 #include "Tpetra_CrsMatrix.hpp" 64 #include "KokkosSparse_CrsMatrix.hpp" 66 #include "Phalanx_DataLayout_MDALayout.hpp" 68 #include "Teuchos_FancyOStream.hpp" 70 template <
typename EvalT,
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
73 const Teuchos::ParameterList& p)
75 std::string scatterName = p.get<std::string>(
"Scatter Name");
76 Teuchos::RCP<PHX::FieldTag> scatterHolder =
77 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
80 const std::vector<std::string>& names =
81 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
83 Teuchos::RCP<PHX::DataLayout> dl =
84 p.get< Teuchos::RCP<const panzer::PureBasis> >(
"Basis")->functional;
87 for (std::size_t eq = 0; eq < names.size(); ++eq) {
88 PHX::MDField<const ScalarT,Cell,NODE>
field = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
91 this->addDependentField(
field.fieldTag());
95 this->addEvaluatedField(*scatterHolder);
97 this->setName(scatterName+
" Scatter Residual");
104 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
107 const Teuchos::ParameterList& p)
108 : globalIndexer_(indexer)
109 , globalDataKey_(
"Residual Scatter Container")
111 std::string scatterName = p.get<std::string>(
"Scatter Name");
113 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
116 const std::vector<std::string>& names =
117 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
120 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
122 Teuchos::RCP<PHX::DataLayout> dl =
123 p.get< Teuchos::RCP<const panzer::PureBasis> >(
"Basis")->functional;
126 scatterFields_.resize(names.size());
127 for (std::size_t eq = 0; eq < names.size(); ++eq) {
128 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
131 this->addDependentField(scatterFields_[eq]);
135 this->addEvaluatedField(*scatterHolder_);
137 if (p.isType<std::string>(
"Global Data Key"))
138 globalDataKey_ = p.get<std::string>(
"Global Data Key");
140 this->setName(scatterName+
" Scatter Residual");
144 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
149 const Workset & workset_0 = (*d.worksets_)[0];
150 const std::string blockId = this->wda(workset_0).
block_id;
152 fieldIds_.resize(scatterFields_.size());
153 fieldOffsets_.resize(scatterFields_.size());
154 fieldGlobalIndexers_.resize(scatterFields_.size());
155 productVectorBlockIndex_.resize(scatterFields_.size());
156 int maxElementBlockGIDCount = -1;
157 for(std::size_t fd=0; fd < scatterFields_.size(); ++fd) {
158 const std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
159 const int globalFieldNum = globalIndexer_->getFieldNum(fieldName);
160 productVectorBlockIndex_[fd] = globalIndexer_->getFieldBlock(globalFieldNum);
161 fieldGlobalIndexers_[fd] = globalIndexer_->getFieldDOFManagers()[productVectorBlockIndex_[fd]];
162 fieldIds_[fd] = fieldGlobalIndexers_[fd]->getFieldNum(fieldName);
164 const std::vector<int>&
offsets = fieldGlobalIndexers_[fd]->getGIDFieldOffsets(blockId,fieldIds_[fd]);
165 fieldOffsets_[fd] = PHX::View<int*>(
"ScatterResidual_BlockedTpetra(Residual):fieldOffsets",
offsets.size());
166 auto hostOffsets = Kokkos::create_mirror_view(fieldOffsets_[fd]);
167 for (std::size_t i=0; i <
offsets.size(); ++i)
169 Kokkos::deep_copy(fieldOffsets_[fd], hostOffsets);
171 maxElementBlockGIDCount = std::max(fieldGlobalIndexers_[fd]->getElementBlockGIDCount(blockId),maxElementBlockGIDCount);
177 worksetLIDs_ = PHX::View<LO**>(
"ScatterResidual_BlockedTpetra(Residual):worksetLIDs",
178 scatterFields_[0].extent(0),
179 maxElementBlockGIDCount);
183 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
188 blockedContainer_ = Teuchos::rcp_dynamic_cast<
const ContainerType>(d.gedc->getDataObject(globalDataKey_),
true);
192 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
197 using Teuchos::rcp_dynamic_cast;
201 const auto& localCellIds = this->wda(workset).cell_local_ids_k;
202 const RCP<ProductVectorBase<double>> thyraBlockResidual = rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer_->get_f(),
true);
205 int currentWorksetLIDSubBlock = -1;
206 for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
208 if (productVectorBlockIndex_[fieldIndex] != currentWorksetLIDSubBlock) {
209 fieldGlobalIndexers_[fieldIndex]->getElementLIDs(localCellIds,worksetLIDs_);
210 currentWorksetLIDSubBlock = productVectorBlockIndex_[fieldIndex];
213 auto& tpetraResidual = *((rcp_dynamic_cast<Thyra::TpetraVector<RealType,LO,GO,NodeT>>(thyraBlockResidual->getNonconstVectorBlock(productVectorBlockIndex_[fieldIndex]),
true))->getTpetraVector());
214 const auto& kokkosResidual = tpetraResidual.getLocalViewDevice(Tpetra::Access::OverwriteAll);
217 const auto& fieldOffsets = fieldOffsets_[fieldIndex];
218 const auto& worksetLIDs = worksetLIDs_;
219 const auto& fieldValues = scatterFields_[fieldIndex].get_static_view();
221 Kokkos::parallel_for(Kokkos::RangePolicy<PHX::Device>(0,workset.num_cells), KOKKOS_LAMBDA (
const int& cell) {
222 for(
int basis=0; basis < static_cast<int>(fieldOffsets.size()); ++basis) {
223 const int lid = worksetLIDs(cell,fieldOffsets(basis));
224 Kokkos::atomic_add(&kokkosResidual(lid,0), fieldValues(cell,basis));
234 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
237 const Teuchos::ParameterList& p)
238 : globalIndexer_(indexer)
239 , globalDataKey_(
"Residual Scatter Container")
241 std::string scatterName = p.get<std::string>(
"Scatter Name");
243 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
246 const std::vector<std::string>& names =
247 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
250 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
252 Teuchos::RCP<PHX::DataLayout> dl =
253 p.get< Teuchos::RCP<const panzer::PureBasis> >(
"Basis")->functional;
256 scatterFields_.resize(names.size());
257 for (std::size_t eq = 0; eq < names.size(); ++eq) {
258 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
261 this->addDependentField(scatterFields_[eq]);
265 this->addEvaluatedField(*scatterHolder_);
267 if (p.isType<std::string>(
"Global Data Key"))
268 globalDataKey_ = p.get<std::string>(
"Global Data Key");
270 this->setName(scatterName+
" Scatter Residual (Jacobian)");
274 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
279 const Workset & workset_0 = (*d.worksets_)[0];
280 const std::string blockId = this->wda(workset_0).
block_id;
282 fieldIds_.resize(scatterFields_.size());
283 fieldOffsets_.resize(scatterFields_.size());
284 productVectorBlockIndex_.resize(scatterFields_.size());
285 for (std::size_t fd=0; fd < scatterFields_.size(); ++fd) {
286 const std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
287 const int globalFieldNum = globalIndexer_->getFieldNum(fieldName);
288 productVectorBlockIndex_[fd] = globalIndexer_->getFieldBlock(globalFieldNum);
289 const auto& fieldGlobalIndexer = globalIndexer_->getFieldDOFManagers()[productVectorBlockIndex_[fd]];
290 fieldIds_[fd] = fieldGlobalIndexer->getFieldNum(fieldName);
292 const std::vector<int>&
offsets = globalIndexer_->getGIDFieldOffsets(blockId,globalFieldNum);
293 fieldOffsets_[fd] = PHX::View<int*>(
"ScatterResidual_BlockedTpetra(Jacobian):fieldOffsets",
offsets.size());
294 auto hostOffsets = Kokkos::create_mirror_view(fieldOffsets_[fd]);
295 for (std::size_t i=0; i <
offsets.size(); ++i)
297 Kokkos::deep_copy(fieldOffsets_[fd], hostOffsets);
303 int elementBlockGIDCount = 0;
304 for (
const auto blockDOFMgr : globalIndexer_->getFieldDOFManagers())
305 elementBlockGIDCount += blockDOFMgr->getElementBlockGIDCount(blockId);
307 worksetLIDs_ = PHX::View<LO**>(
"ScatterResidual_BlockedTpetra(Jacobian):worksetLIDs",
308 scatterFields_[0].extent(0),
309 elementBlockGIDCount);
312 const auto& blockGlobalIndexers = globalIndexer_->getFieldDOFManagers();
313 const int numBlocks =
static_cast<int>(globalIndexer_->getFieldDOFManagers().size());
314 blockOffsets_ = PHX::View<LO*>(
"ScatterResidual_BlockedTpetra(Jacobian):blockOffsets_",
316 const auto hostBlockOffsets = Kokkos::create_mirror_view(blockOffsets_);
317 for (
int blk=0;blk<numBlocks;blk++) {
318 int blockOffset = globalIndexer_->getBlockGIDOffset(blockId,blk);
319 hostBlockOffsets(blk) = blockOffset;
321 hostBlockOffsets(numBlocks) = hostBlockOffsets(numBlocks-1) + blockGlobalIndexers[blockGlobalIndexers.size()-1]->getElementBlockGIDCount(blockId);
322 Kokkos::deep_copy(blockOffsets_,hostBlockOffsets);
327 for (
int blk=0;blk<numBlocks;blk++) {
328 const int blockDerivativeSize = hostBlockOffsets(blk+1) - hostBlockOffsets(blk);
329 TEUCHOS_TEST_FOR_EXCEPTION(blockDerivativeSize > 256, std::runtime_error,
330 "ERROR: the derivative dimension for sub block " 331 << blk <<
"with a value of " << blockDerivativeSize
332 <<
"is larger than the size allocated for cLIDs and vals " 333 <<
"in the evaluate call! You must manually increase the " 334 <<
"size and recompile!");
339 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
344 using Teuchos::rcp_dynamic_cast;
347 blockedContainer_ = rcp_dynamic_cast<
const ContainerType>(d.gedc->getDataObject(globalDataKey_));
349 if(blockedContainer_==Teuchos::null) {
350 RCP<const LOCPair_GlobalEvaluationData> gdata = rcp_dynamic_cast<
const LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),
true);
351 blockedContainer_ = rcp_dynamic_cast<
const ContainerType>(gdata->getGhostedLOC());
356 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
361 using Teuchos::rcp_dynamic_cast;
366 const auto& localCellIds = this->wda(workset).cell_local_ids_k;
368 const int numFieldBlocks = globalIndexer_->getNumFieldBlocks();
369 const RCP<const ContainerType> blockedContainer = blockedContainer_;
370 const RCP<ProductVectorBase<double>> thyraBlockResidual = rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer_->get_f());
371 const bool haveResidual = Teuchos::nonnull(thyraBlockResidual);
372 const RCP<BlockedLinearOpBase<double>> Jac = rcp_dynamic_cast<BlockedLinearOpBase<double> >(blockedContainer_->get_A(),
true);
378 using LocalMatrixType = KokkosSparse::CrsMatrix<double,LO,PHX::Device,Kokkos::MemoryTraits<Kokkos::Unmanaged>,
size_t>;
379 typename PHX::View<LocalMatrixType**>::HostMirror
380 hostJacTpetraBlocks(
"panzer::ScatterResidual_BlockTpetra<Jacobian>::hostJacTpetraBlocks", numFieldBlocks,numFieldBlocks);
382 PHX::View<int**> blockExistsInJac = PHX::View<int**>(
"blockExistsInJac_",numFieldBlocks,numFieldBlocks);
383 auto hostBlockExistsInJac = Kokkos::create_mirror_view(blockExistsInJac);
385 for (
int row=0; row < numFieldBlocks; ++row) {
386 for (
int col=0; col < numFieldBlocks; ++col) {
387 const auto thyraTpetraOperator = rcp_dynamic_cast<Thyra::TpetraLinearOp<double,LO,GO,NodeT>>(Jac->getNonconstBlock(row,col),
false);
388 if (nonnull(thyraTpetraOperator)) {
404 const auto tpetraCrsMatrix = rcp_dynamic_cast<Tpetra::CrsMatrix<double,LO,GO,NodeT>>(thyraTpetraOperator->getTpetraOperator(),
true);
405 const auto managedMatrix = tpetraCrsMatrix->getLocalMatrixDevice();
406 const auto managedGraph = managedMatrix.graph;
409 using StaticCrsGraphType =
typename LocalMatrixType::StaticCrsGraphType;
410 StaticCrsGraphType unmanagedGraph;
411 unmanagedGraph.entries =
typename StaticCrsGraphType::entries_type(managedGraph.entries.data(),managedGraph.entries.extent(0));
412 unmanagedGraph.row_map =
typename StaticCrsGraphType::row_map_type(managedGraph.row_map.data(),managedGraph.row_map.extent(0));
413 unmanagedGraph.row_block_offsets =
typename StaticCrsGraphType::row_block_type(managedGraph.row_block_offsets.data(),managedGraph.row_block_offsets.extent(0));
415 typename LocalMatrixType::values_type unmanagedValues(managedMatrix.values.data(),managedMatrix.values.extent(0));
416 LocalMatrixType unmanagedMatrix(managedMatrix.values.label(), managedMatrix.numCols(), unmanagedValues, unmanagedGraph);
417 new (&hostJacTpetraBlocks(row,col)) LocalMatrixType(unmanagedMatrix);
420 hostBlockExistsInJac(row,col) = 1;
423 hostBlockExistsInJac(row,col) = 0;
427 typename PHX::View<LocalMatrixType**>
428 jacTpetraBlocks(
"panzer::ScatterResidual_BlockedTpetra<Jacobian>::jacTpetraBlocks",numFieldBlocks,numFieldBlocks);
429 Kokkos::deep_copy(jacTpetraBlocks,hostJacTpetraBlocks);
430 Kokkos::deep_copy(blockExistsInJac,hostBlockExistsInJac);
437 const auto& globalIndexers = globalIndexer_->getFieldDOFManagers();
438 auto blockOffsets_h = Kokkos::create_mirror_view(blockOffsets_);
439 Kokkos::deep_copy(blockOffsets_h, blockOffsets_);
440 for (
size_t block=0; block < globalIndexers.size(); ++block) {
441 const auto subviewOfBlockLIDs = Kokkos::subview(worksetLIDs_,Kokkos::ALL(), std::make_pair(blockOffsets_h(block),blockOffsets_h(block+1)));
442 globalIndexers[block]->getElementLIDs(localCellIds,subviewOfBlockLIDs);
446 for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
448 const int blockRowIndex = productVectorBlockIndex_[fieldIndex];
449 typename Tpetra::Vector<double,LO,GO,PHX::Device>::dual_view_type::t_dev kokkosResidual;
451 auto& tpetraResidual = *((rcp_dynamic_cast<Thyra::TpetraVector<RealType,LO,GO,NodeT>>(thyraBlockResidual->getNonconstVectorBlock(blockRowIndex),
true))->getTpetraVector());
452 kokkosResidual = tpetraResidual.getLocalViewDevice(Tpetra::Access::OverwriteAll);
456 const PHX::View<const int*> fieldOffsets = fieldOffsets_[fieldIndex];
457 const PHX::View<const LO**> worksetLIDs = worksetLIDs_;
458 const PHX::View<const ScalarT**> fieldValues = scatterFields_[fieldIndex].get_static_view();
459 const PHX::View<const LO*> blockOffsets = blockOffsets_;
461 Kokkos::parallel_for(Kokkos::RangePolicy<PHX::Device>(0,workset.num_cells), KOKKOS_LAMBDA (
const int& cell) {
463 typename Sacado::ScalarType<ScalarT>::type vals[256];
465 for(
int basis=0; basis < static_cast<int>(fieldOffsets.size()); ++basis) {
466 typedef PHX::MDField<const ScalarT,Cell,NODE>
FieldType;
467 typename FieldType::array_type::reference_type tmpFieldVal = fieldValues(cell,basis);
468 const int rowLID = worksetLIDs(cell,fieldOffsets(basis));
471 Kokkos::atomic_add(&kokkosResidual(rowLID,0), tmpFieldVal.val());
473 for (
int blockColIndex=0; blockColIndex < numFieldBlocks; ++blockColIndex) {
474 if (blockExistsInJac(blockRowIndex,blockColIndex)) {
475 const int start = blockOffsets(blockColIndex);
476 const int stop = blockOffsets(blockColIndex+1);
477 const int sensSize = stop-start;
479 for (
int i=0; i < sensSize; ++i) {
480 cLIDs[i] = worksetLIDs(cell,start+i);
481 vals[i] = tmpFieldVal.fastAccessDx(start+i);
483 jacTpetraBlocks(blockRowIndex,blockColIndex).sumIntoValues(rowLID,cLIDs,sensSize,vals,
true,
true);
492 for (
int row=0; row < numFieldBlocks; ++row) {
493 for (
int col=0; col < numFieldBlocks; ++col) {
494 if (hostBlockExistsInJac(row,col)) {
495 hostJacTpetraBlocks(row,col).~CrsMatrix();
Pushes residual values into the residual vector for a Newton-based solve.
void postRegistrationSetup(typename TRAITS::SetupData, PHX::FieldManager< TRAITS > &)
ScatterResidual_BlockedTpetra(const Teuchos::RCP< const BlockedDOFManager > &)
FieldType
The type of discretization to use for a field pattern.
PHX::View< const int * > offsets
void evaluateFields(typename TRAITS::EvalData)
std::string block_id
DEPRECATED - use: getElementBlock()
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...