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"
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>
75 std::string scatterName = p.
get<std::string>(
"Scatter Name");
80 const std::vector<std::string>& names =
87 for (std::size_t eq = 0; eq < names.size(); ++eq) {
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>
108 : globalIndexer_(indexer)
109 , globalDataKey_(
"Residual Scatter Container")
111 std::string scatterName = p.
get<std::string>(
"Scatter Name");
116 const std::vector<std::string>& names =
126 scatterFields_.resize(names.size());
127 for (std::size_t eq = 0; eq < names.size(); ++eq) {
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)
168 hostOffsets(i) = offsets[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 using Teuchos::rcp_dynamic_cast;
191 blockedContainer_ = rcp_dynamic_cast<
const ContainerType>(d.gedc->getDataObject(globalDataKey_));
193 if(blockedContainer_==Teuchos::null) {
200 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
205 using Teuchos::rcp_dynamic_cast;
209 const auto& localCellIds = this->wda(workset).cell_local_ids_k;
213 int currentWorksetLIDSubBlock = -1;
214 for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
216 if (productVectorBlockIndex_[fieldIndex] != currentWorksetLIDSubBlock) {
217 fieldGlobalIndexers_[fieldIndex]->getElementLIDs(localCellIds,worksetLIDs_);
218 currentWorksetLIDSubBlock = productVectorBlockIndex_[fieldIndex];
221 auto& tpetraResidual = *((rcp_dynamic_cast<Thyra::TpetraVector<RealType,LO,GO,NodeT>>(thyraBlockResidual->getNonconstVectorBlock(productVectorBlockIndex_[fieldIndex]),
true))->getTpetraVector());
222 const auto& kokkosResidual = tpetraResidual.getLocalViewDevice(Tpetra::Access::ReadWrite);
225 const auto& fieldOffsets = fieldOffsets_[fieldIndex];
226 const auto& worksetLIDs = worksetLIDs_;
227 const auto& fieldValues = scatterFields_[fieldIndex].get_static_view();
229 Kokkos::parallel_for(Kokkos::RangePolicy<PHX::Device>(0,workset.num_cells), KOKKOS_LAMBDA (
const int& cell) {
230 for(
int basis=0; basis < static_cast<int>(fieldOffsets.size()); ++basis) {
231 const int lid = worksetLIDs(cell,fieldOffsets(basis));
232 Kokkos::atomic_add(&kokkosResidual(lid,0), fieldValues(cell,basis));
242 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
246 : globalIndexer_(indexer)
247 , globalDataKey_(
"Residual Scatter Container")
249 std::string scatterName = p.
get<std::string>(
"Scatter Name");
254 const std::vector<std::string>& names =
264 scatterFields_.resize(names.size());
265 for (std::size_t eq = 0; eq < names.size(); ++eq) {
269 this->addDependentField(scatterFields_[eq]);
273 this->addEvaluatedField(*scatterHolder_);
275 if (p.
isType<std::string>(
"Global Data Key"))
276 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
278 this->setName(scatterName+
" Scatter Residual (Jacobian)");
282 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
287 const Workset & workset_0 = (*d.worksets_)[0];
288 const std::string blockId = this->wda(workset_0).
block_id;
290 fieldIds_.resize(scatterFields_.size());
291 fieldOffsets_.resize(scatterFields_.size());
292 productVectorBlockIndex_.resize(scatterFields_.size());
293 for (std::size_t fd=0; fd < scatterFields_.size(); ++fd) {
294 const std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
295 const int globalFieldNum = globalIndexer_->getFieldNum(fieldName);
296 productVectorBlockIndex_[fd] = globalIndexer_->getFieldBlock(globalFieldNum);
297 const auto& fieldGlobalIndexer = globalIndexer_->getFieldDOFManagers()[productVectorBlockIndex_[fd]];
298 fieldIds_[fd] = fieldGlobalIndexer->getFieldNum(fieldName);
300 const std::vector<int>&
offsets = globalIndexer_->getGIDFieldOffsets(blockId,globalFieldNum);
301 fieldOffsets_[fd] = PHX::View<int*>(
"ScatterResidual_BlockedTpetra(Jacobian):fieldOffsets",offsets.size());
302 auto hostOffsets = Kokkos::create_mirror_view(fieldOffsets_[fd]);
303 for (std::size_t i=0; i < offsets.size(); ++i)
304 hostOffsets(i) = offsets[i];
305 Kokkos::deep_copy(fieldOffsets_[fd], hostOffsets);
311 int elementBlockGIDCount = 0;
312 for (
const auto& blockDOFMgr : globalIndexer_->getFieldDOFManagers())
313 elementBlockGIDCount += blockDOFMgr->getElementBlockGIDCount(blockId);
316 "ScatterResidual_BlockedTpetra(Jacobian):worksetLIDs_",
317 scatterFields_[0].extent(0), elementBlockGIDCount );
320 const auto& blockGlobalIndexers = globalIndexer_->getFieldDOFManagers();
321 const int numBlocks =
static_cast<int>(globalIndexer_->getFieldDOFManagers().size());
322 blockOffsets_ = PHX::View<LO*>(
"ScatterResidual_BlockedTpetra(Jacobian):blockOffsets_",
324 const auto hostBlockOffsets = Kokkos::create_mirror_view(blockOffsets_);
325 for (
int blk=0;blk<numBlocks;blk++) {
326 int blockOffset = globalIndexer_->getBlockGIDOffset(blockId,blk);
327 hostBlockOffsets(blk) = blockOffset;
329 hostBlockOffsets(numBlocks) = hostBlockOffsets(numBlocks-1) + blockGlobalIndexers[blockGlobalIndexers.size()-1]->getElementBlockGIDCount(blockId);
330 Kokkos::deep_copy(blockOffsets_,hostBlockOffsets);
334 int max_blockDerivativeSize = 0;
335 for (
int blk=0;blk<numBlocks;blk++) {
336 const int blockDerivativeSize = hostBlockOffsets(blk+1) - hostBlockOffsets(blk);
337 if ( blockDerivativeSize > max_blockDerivativeSize )
338 max_blockDerivativeSize = blockDerivativeSize;
341 "ScatterResidual_BlockedTpetra(Jacobian):workset_vals_",
342 scatterFields_[0].extent(0), max_blockDerivativeSize );
346 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
351 using Teuchos::rcp_dynamic_cast;
354 blockedContainer_ = rcp_dynamic_cast<
const ContainerType>(d.gedc->getDataObject(globalDataKey_));
356 if(blockedContainer_==Teuchos::null) {
363 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
368 using Teuchos::rcp_dynamic_cast;
373 const auto& localCellIds = this->wda(workset).cell_local_ids_k;
375 const int numFieldBlocks = globalIndexer_->getNumFieldBlocks();
385 using LocalMatrixType = KokkosSparse::CrsMatrix<double,LO,PHX::Device,Kokkos::MemoryTraits<Kokkos::Unmanaged>,
size_t>;
386 typename PHX::View<LocalMatrixType**>::HostMirror
387 hostJacTpetraBlocks(
"panzer::ScatterResidual_BlockTpetra<Jacobian>::hostJacTpetraBlocks", numFieldBlocks,numFieldBlocks);
389 PHX::View<int**> blockExistsInJac = PHX::View<int**>(
"blockExistsInJac_",numFieldBlocks,numFieldBlocks);
390 auto hostBlockExistsInJac = Kokkos::create_mirror_view(blockExistsInJac);
392 for (
int row=0; row < numFieldBlocks; ++row) {
393 for (
int col=0; col < numFieldBlocks; ++col) {
394 const auto thyraTpetraOperator = rcp_dynamic_cast<Thyra::TpetraLinearOp<double,LO,GO,NodeT>>(Jac->getNonconstBlock(row,col),
false);
395 if (
nonnull(thyraTpetraOperator)) {
411 const auto tpetraCrsMatrix = rcp_dynamic_cast<Tpetra::CrsMatrix<double,LO,GO,NodeT>>(thyraTpetraOperator->getTpetraOperator(),
true);
412 const auto managedMatrix = tpetraCrsMatrix->getLocalMatrixDevice();
413 const auto managedGraph = managedMatrix.graph;
416 using StaticCrsGraphType =
typename LocalMatrixType::StaticCrsGraphType;
417 StaticCrsGraphType unmanagedGraph;
418 unmanagedGraph.entries =
typename StaticCrsGraphType::entries_type(managedGraph.entries.data(),managedGraph.entries.extent(0));
419 unmanagedGraph.row_map =
typename StaticCrsGraphType::row_map_type(managedGraph.row_map.data(),managedGraph.row_map.extent(0));
420 unmanagedGraph.row_block_offsets =
typename StaticCrsGraphType::row_block_type(managedGraph.row_block_offsets.data(),managedGraph.row_block_offsets.extent(0));
422 typename LocalMatrixType::values_type unmanagedValues(managedMatrix.values.data(),managedMatrix.values.extent(0));
423 LocalMatrixType unmanagedMatrix(managedMatrix.values.label(), managedMatrix.numCols(), unmanagedValues, unmanagedGraph);
424 new (&hostJacTpetraBlocks(row,col)) LocalMatrixType(unmanagedMatrix);
427 hostBlockExistsInJac(row,col) = 1;
430 hostBlockExistsInJac(row,col) = 0;
434 typename PHX::View<LocalMatrixType**>
435 jacTpetraBlocks(
"panzer::ScatterResidual_BlockedTpetra<Jacobian>::jacTpetraBlocks",numFieldBlocks,numFieldBlocks);
436 Kokkos::deep_copy(jacTpetraBlocks,hostJacTpetraBlocks);
437 Kokkos::deep_copy(blockExistsInJac,hostBlockExistsInJac);
444 const auto& globalIndexers = globalIndexer_->getFieldDOFManagers();
445 auto blockOffsets_h = Kokkos::create_mirror_view(blockOffsets_);
446 Kokkos::deep_copy(blockOffsets_h, blockOffsets_);
447 for (
size_t block=0; block < globalIndexers.size(); ++block) {
448 const auto subviewOfBlockLIDs = Kokkos::subview(worksetLIDs_,Kokkos::ALL(), std::make_pair(blockOffsets_h(block),blockOffsets_h(block+1)));
449 globalIndexers[block]->getElementLIDs(localCellIds,subviewOfBlockLIDs);
453 for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
455 const int blockRowIndex = productVectorBlockIndex_[fieldIndex];
456 typename Tpetra::Vector<double,LO,GO,PHX::Device>::dual_view_type::t_dev kokkosResidual;
458 auto& tpetraResidual = *((rcp_dynamic_cast<Thyra::TpetraVector<RealType,LO,GO,NodeT>>(thyraBlockResidual->getNonconstVectorBlock(blockRowIndex),
true))->getTpetraVector());
459 kokkosResidual = tpetraResidual.getLocalViewDevice(Tpetra::Access::ReadWrite);
463 const PHX::View<const int*> fieldOffsets = fieldOffsets_[fieldIndex];
464 const PHX::View<const ScalarT**> fieldValues = scatterFields_[fieldIndex].get_static_view();
465 const PHX::View<const LO*> blockOffsets = blockOffsets_;
470 Kokkos::parallel_for(Kokkos::RangePolicy<PHX::Device>(0,workset.num_cells), KOKKOS_LAMBDA (
const int& cell) {
471 for(
int basis=0; basis < static_cast<int>(fieldOffsets.size()); ++basis) {
473 typename FieldType::array_type::reference_type tmpFieldVal = fieldValues(cell,basis);
474 const int rowLID = worksetLIDs(cell,fieldOffsets(basis));
477 Kokkos::atomic_add(&kokkosResidual(rowLID,0), tmpFieldVal.val());
479 for (
int blockColIndex=0; blockColIndex < numFieldBlocks; ++blockColIndex) {
480 if (blockExistsInJac(blockRowIndex,blockColIndex)) {
481 const int start = blockOffsets(blockColIndex);
482 const int stop = blockOffsets(blockColIndex+1);
483 const int sensSize = stop-start;
485 for (
int i=0; i < sensSize; ++i)
486 workset_vals(cell,i) = tmpFieldVal.fastAccessDx(start+i);
488 jacTpetraBlocks(blockRowIndex,blockColIndex).sumIntoValues(rowLID, &worksetLIDs(cell,start), sensSize, &workset_vals(cell,0),
true,
true);
497 for (
int row=0; row < numFieldBlocks; ++row) {
498 for (
int col=0; col < numFieldBlocks; ++col) {
499 if (hostBlockExistsInJac(row,col)) {
500 hostJacTpetraBlocks(row,col).~CrsMatrix();
T & get(const std::string &name, T def_value)
bool nonnull(const std::shared_ptr< T > &p)
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.
Teuchos::RCP< LinearObjContainer > getGhostedLOC() const
PHX::View< const int * > offsets
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void evaluateFields(typename TRAITS::EvalData)
virtual void preEvaluate(typename Traits::PreEvalData d)=0
bool nonnull(const boost::shared_ptr< T > &p)
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...
bool isType(const std::string &name) const