43 #ifndef PANZER_SCATTER_DIRICHLET_RESIDUAL_BLOCEDTPETRA_IMPL_HPP
44 #define PANZER_SCATTER_DIRICHLET_RESIDUAL_BLOCEDTPETRA_IMPL_HPP
47 #include "Teuchos_Assert.hpp"
49 #include "Phalanx_DataLayout.hpp"
61 #include "Phalanx_DataLayout_MDALayout.hpp"
63 #include "Thyra_SpmdVectorBase.hpp"
64 #include "Thyra_ProductVectorBase.hpp"
65 #include "Thyra_BlockedLinearOpBase.hpp"
68 #include "Teuchos_FancyOStream.hpp"
70 #include <unordered_map>
72 template <
typename EvalT,
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
77 std::string scatterName = p.
get<std::string>(
"Scatter Name");
82 const std::vector<std::string>& names =
89 for (std::size_t eq = 0; eq < names.size(); ++eq) {
93 this->addDependentField(scatterField.fieldTag());
97 this->addEvaluatedField(*scatterHolder);
99 this->setName(scatterName+
" Scatter Residual");
107 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
111 : globalIndexer_(indexer)
112 , globalDataKey_(
"Residual Scatter Container")
114 std::string scatterName = p.
get<std::string>(
"Scatter Name");
119 const std::vector<std::string>& names =
126 scatterIC_ = p.
isParameter(
"Scatter Initial Condition") ? p.
get<
bool>(
"Scatter Initial Condition") :
false;
132 side_subcell_dim_ = p.
get<
int>(
"Side Subcell Dimension");
133 local_side_id_ = p.
get<
int>(
"Local Side ID");
137 scatterFields_.resize(names.size());
138 for (std::size_t eq = 0; eq < names.size(); ++eq) {
142 this->addDependentField(scatterFields_[eq]);
145 checkApplyBC_ = p.
isParameter(
"Check Apply BC") ? p.
get<
bool>(
"Check Apply BC") :
false;
146 applyBC_.resize(names.size());
148 for (std::size_t eq = 0; eq < names.size(); ++eq) {
150 this->addDependentField(applyBC_[eq]);
155 this->addEvaluatedField(*scatterHolder_);
157 if (p.
isType<std::string>(
"Global Data Key"))
158 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
160 this->setName(scatterName+
" Scatter Residual");
164 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
169 const Workset & workset_0 = (*d.worksets_)[0];
170 const std::string blockId = this->wda(workset_0).
block_id;
172 fieldIds_.resize(scatterFields_.size());
173 fieldOffsets_.resize(scatterFields_.size());
174 basisIndexForMDFieldOffsets_.resize(scatterFields_.size());
175 fieldGlobalIndexers_.resize(scatterFields_.size());
176 productVectorBlockIndex_.resize(scatterFields_.size());
177 int maxElementBlockGIDCount = -1;
178 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
180 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
182 const int globalFieldNum = globalIndexer_->getFieldNum(fieldName);
183 productVectorBlockIndex_[fd] = globalIndexer_->getFieldBlock(globalFieldNum);
184 fieldGlobalIndexers_[fd] = globalIndexer_->getFieldDOFManagers()[productVectorBlockIndex_[fd]];
185 fieldIds_[fd] = fieldGlobalIndexers_[fd]->getFieldNum(fieldName);
189 const auto& offsetPair = fieldGlobalIndexers_[fd]->getGIDFieldOffsets_closure(blockId,fieldIds_[fd],side_subcell_dim_,local_side_id_);
191 const auto&
offsets = offsetPair.first;
192 fieldOffsets_[fd] = PHX::View<int*>(
"ScatterDirichletResidual_BlockedTpetra(Residual):fieldOffsets",
offsets.size());
193 auto hostOffsets = Kokkos::create_mirror_view(fieldOffsets_[fd]);
194 for (std::size_t i=0; i <
offsets.size(); ++i)
196 Kokkos::deep_copy(fieldOffsets_[fd], hostOffsets);
199 const auto& basisIndex = offsetPair.second;
200 basisIndexForMDFieldOffsets_[fd] = PHX::View<int*>(
"ScatterDirichletResidual_BlockedTpetra(Residual):basisIndexForMDFieldOffsets",basisIndex.size());
201 auto hostBasisIndex = Kokkos::create_mirror_view(basisIndexForMDFieldOffsets_[fd]);
202 for (std::size_t i=0; i < basisIndex.size(); ++i)
203 hostBasisIndex(i) = basisIndex[i];
204 Kokkos::deep_copy(basisIndexForMDFieldOffsets_[fd], hostBasisIndex);
209 const std::vector<int>&
offsets = fieldGlobalIndexers_[fd]->getGIDFieldOffsets(blockId,fieldIds_[fd]);
210 fieldOffsets_[fd] = PHX::View<int*>(
"ScatterDirichletResidual_BlockedTpetra(Residual):fieldOffsets",offsets.size());
211 auto hostOffsets = Kokkos::create_mirror_view(fieldOffsets_[fd]);
212 for (std::size_t i=0; i < offsets.size(); ++i)
213 hostOffsets(i) = offsets[i];
214 Kokkos::deep_copy(fieldOffsets_[fd], hostOffsets);
217 maxElementBlockGIDCount = std::max(fieldGlobalIndexers_[fd]->getElementBlockGIDCount(blockId),maxElementBlockGIDCount);
223 worksetLIDs_ = PHX::View<LO**>(
"ScatterResidual_BlockedTpetra(Residual):worksetLIDs",
224 scatterFields_[0].extent(0),
225 maxElementBlockGIDCount);
229 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
235 = Teuchos::rcp_dynamic_cast<
ContainerType>(d.gedc->getDataObject(
"Dirichlet Counter"),
true);
241 blockedContainer_ = Teuchos::rcp_dynamic_cast<
const ContainerType>(d.gedc->getDataObject(globalDataKey_),
true);
246 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
251 using Teuchos::rcp_dynamic_cast;
255 const auto& localCellIds = this->wda(workset).cell_local_ids_k;
258 rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer_->get_f(),
true) :
259 rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer_->get_x(),
true);
262 int currentWorksetLIDSubBlock = -1;
263 for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
265 if (productVectorBlockIndex_[fieldIndex] != currentWorksetLIDSubBlock) {
266 fieldGlobalIndexers_[fieldIndex]->getElementLIDs(localCellIds,worksetLIDs_);
267 currentWorksetLIDSubBlock = productVectorBlockIndex_[fieldIndex];
271 auto& tpetraScatterTarget = *((rcp_dynamic_cast<Thyra::TpetraVector<RealType,LO,GO,NodeT>>(thyraScatterTarget->getNonconstVectorBlock(productVectorBlockIndex_[fieldIndex]),
true))->getTpetraVector());
272 const auto& kokkosScatterTarget = tpetraScatterTarget.getLocalViewDevice(Tpetra::Access::ReadWrite);
275 auto& tpetraDirichletCounter = *((rcp_dynamic_cast<Thyra::TpetraVector<RealType,LO,GO,NodeT>>(dirichletCounter_->getNonconstVectorBlock(productVectorBlockIndex_[fieldIndex]),
true))->getTpetraVector());
276 const auto& kokkosDirichletCounter = tpetraDirichletCounter.getLocalViewDevice(Tpetra::Access::ReadWrite);
279 const auto fieldOffsets = fieldOffsets_[fieldIndex];
280 const auto basisIndices = basisIndexForMDFieldOffsets_[fieldIndex];
281 const auto worksetLIDs = worksetLIDs_;
282 const auto fieldValues = scatterFields_[fieldIndex].get_static_view();
283 const auto applyBC = applyBC_[fieldIndex].get_static_view();
284 const bool checkApplyBC = checkApplyBC_;
288 Kokkos::parallel_for(Kokkos::RangePolicy<PHX::Device>(0,workset.num_cells), KOKKOS_LAMBDA (
const int& cell) {
289 for (
int basis=0; basis < static_cast<int>(fieldOffsets.size()); ++basis) {
290 const int lid = worksetLIDs(cell,fieldOffsets(basis));
293 const int basisIndex = basisIndices(basis);
297 if (!applyBC(cell,basisIndex))
300 kokkosScatterTarget(lid,0) = fieldValues(cell,basisIndex);
301 kokkosDirichletCounter(lid,0) = 1.0;
307 Kokkos::parallel_for(Kokkos::RangePolicy<PHX::Device>(0,workset.num_cells), KOKKOS_LAMBDA (
const int& cell) {
308 for (
int basis=0; basis < static_cast<int>(fieldOffsets.size()); ++basis) {
309 const int lid = worksetLIDs(cell,fieldOffsets(basis));
312 kokkosScatterTarget(lid,0) = fieldValues(cell,basis);
313 kokkosDirichletCounter(lid,0) = 1.0;
326 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
330 : globalIndexer_(indexer)
331 , globalDataKey_(
"Residual Scatter Container")
333 std::string scatterName = p.
get<std::string>(
"Scatter Name");
338 const std::vector<std::string>& names =
347 side_subcell_dim_ = p.
get<
int>(
"Side Subcell Dimension");
348 local_side_id_ = p.
get<
int>(
"Local Side ID");
351 scatterFields_.resize(names.size());
352 for (std::size_t eq = 0; eq < names.size(); ++eq) {
356 this->addDependentField(scatterFields_[eq]);
359 checkApplyBC_ = p.
get<
bool>(
"Check Apply BC");
360 applyBC_.resize(names.size());
362 for (std::size_t eq = 0; eq < names.size(); ++eq) {
364 this->addDependentField(applyBC_[eq]);
369 this->addEvaluatedField(*scatterHolder_);
371 if (p.
isType<std::string>(
"Global Data Key"))
372 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
374 this->setName(scatterName+
" Scatter Residual (Jacobian)");
378 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
383 const Workset & workset_0 = (*d.worksets_)[0];
384 const std::string blockId = this->wda(workset_0).
block_id;
386 fieldIds_.resize(scatterFields_.size());
387 fieldOffsets_.resize(scatterFields_.size());
388 basisIndexForMDFieldOffsets_.resize(scatterFields_.size());
389 productVectorBlockIndex_.resize(scatterFields_.size());
390 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
392 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
393 const int globalFieldNum = globalIndexer_->getFieldNum(fieldName);
394 productVectorBlockIndex_[fd] = globalIndexer_->getFieldBlock(globalFieldNum);
395 const auto& fieldGlobalIndexer = globalIndexer_->getFieldDOFManagers()[productVectorBlockIndex_[fd]];
396 fieldIds_[fd] = fieldGlobalIndexer->getFieldNum(fieldName);
398 const auto& offsetPair = fieldGlobalIndexer->getGIDFieldOffsets_closure(blockId,fieldIds_[fd],side_subcell_dim_,local_side_id_);
400 const auto&
offsets = offsetPair.first;
401 fieldOffsets_[fd] = PHX::View<int*>(
"ScatterDirichletResidual_BlockedTpetra(Jacobian):fieldOffsets",
offsets.size());
402 auto hostOffsets = Kokkos::create_mirror_view(fieldOffsets_[fd]);
403 for (std::size_t i=0; i <
offsets.size(); ++i)
405 Kokkos::deep_copy(fieldOffsets_[fd], hostOffsets);
408 const auto& basisIndex = offsetPair.second;
409 basisIndexForMDFieldOffsets_[fd] = PHX::View<int*>(
"ScatterDirichletResidual_BlockedTpetra(Jacobian):basisIndexForMDFieldOffsets",basisIndex.size());
410 auto hostBasisIndex = Kokkos::create_mirror_view(basisIndexForMDFieldOffsets_[fd]);
411 for (std::size_t i=0; i < basisIndex.size(); ++i)
412 hostBasisIndex(i) = basisIndex[i];
413 Kokkos::deep_copy(basisIndexForMDFieldOffsets_[fd], hostBasisIndex);
420 int elementBlockGIDCount = 0;
421 for (
const auto& blockDOFMgr : globalIndexer_->getFieldDOFManagers())
422 elementBlockGIDCount += blockDOFMgr->getElementBlockGIDCount(blockId);
424 worksetLIDs_ = PHX::View<LO**>(
"ScatterDirichletResidual_BlockedTpetra(Jacobian):worksetLIDs",
425 scatterFields_[0].extent(0),
426 elementBlockGIDCount);
429 const auto& blockGlobalIndexers = globalIndexer_->getFieldDOFManagers();
430 const int numBlocks =
static_cast<int>(globalIndexer_->getFieldDOFManagers().size());
431 blockOffsets_ = PHX::View<LO*>(
"ScatterDirichletResidual_BlockedTpetra(Jacobian):blockOffsets_",
433 const auto hostBlockOffsets = Kokkos::create_mirror_view(blockOffsets_);
434 for (
int blk=0;blk<numBlocks;blk++) {
435 int blockOffset = globalIndexer_->getBlockGIDOffset(blockId,blk);
436 hostBlockOffsets(blk) = blockOffset;
438 hostBlockOffsets(numBlocks) = hostBlockOffsets(numBlocks-1) + blockGlobalIndexers[blockGlobalIndexers.size()-1]->getElementBlockGIDCount(blockId);
439 Kokkos::deep_copy(blockOffsets_,hostBlockOffsets);
444 for (
int blk=0;blk<numBlocks;blk++) {
445 const int blockDerivativeSize = hostBlockOffsets(blk+1) - hostBlockOffsets(blk);
447 "ERROR: the derivative dimension for sub block "
448 << blk <<
"with a value of " << blockDerivativeSize
449 <<
"is larger than the size allocated for cLIDs and vals "
450 <<
"in the evaluate call! You must manually increase the "
451 <<
"size and recompile!");
456 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
462 = Teuchos::rcp_dynamic_cast<
const ContainerType>(d.gedc->getDataObject(
"Dirichlet Counter"),
true);
468 blockedContainer_ = Teuchos::rcp_dynamic_cast<
const ContainerType>(d.gedc->getDataObject(globalDataKey_),
true);
473 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
478 using Teuchos::rcp_dynamic_cast;
483 const auto& localCellIds = this->wda(workset).cell_local_ids_k;
485 const int numFieldBlocks = globalIndexer_->getNumFieldBlocks();
495 using LocalMatrixType = KokkosSparse::CrsMatrix<double,LO,PHX::Device,Kokkos::MemoryTraits<Kokkos::Unmanaged>,
size_t>;
496 typename PHX::View<LocalMatrixType**>::HostMirror
497 hostJacTpetraBlocks(
"panzer::ScatterResidual_BlockTpetra<Jacobian>::hostJacTpetraBlocks", numFieldBlocks,numFieldBlocks);
499 PHX::View<int**> blockExistsInJac = PHX::View<int**>(
"blockExistsInJac_",numFieldBlocks,numFieldBlocks);
500 auto hostBlockExistsInJac = Kokkos::create_mirror_view(blockExistsInJac);
502 for (
int row=0; row < numFieldBlocks; ++row) {
503 for (
int col=0; col < numFieldBlocks; ++col) {
504 const auto thyraTpetraOperator = rcp_dynamic_cast<Thyra::TpetraLinearOp<double,LO,GO,NodeT>>(Jac->getNonconstBlock(row,col),
false);
505 if (
nonnull(thyraTpetraOperator)) {
521 const auto tpetraCrsMatrix = rcp_dynamic_cast<Tpetra::CrsMatrix<double,LO,GO,NodeT>>(thyraTpetraOperator->getTpetraOperator(),
true);
522 const auto managedMatrix = tpetraCrsMatrix->getLocalMatrixDevice();
523 const auto managedGraph = managedMatrix.graph;
526 using StaticCrsGraphType =
typename LocalMatrixType::StaticCrsGraphType;
527 StaticCrsGraphType unmanagedGraph;
528 unmanagedGraph.entries =
typename StaticCrsGraphType::entries_type(managedGraph.entries.data(),managedGraph.entries.extent(0));
529 unmanagedGraph.row_map =
typename StaticCrsGraphType::row_map_type(managedGraph.row_map.data(),managedGraph.row_map.extent(0));
530 unmanagedGraph.row_block_offsets =
typename StaticCrsGraphType::row_block_type(managedGraph.row_block_offsets.data(),managedGraph.row_block_offsets.extent(0));
532 typename LocalMatrixType::values_type unmanagedValues(managedMatrix.values.data(),managedMatrix.values.extent(0));
533 LocalMatrixType unmanagedMatrix(managedMatrix.values.label(), managedMatrix.numCols(), unmanagedValues, unmanagedGraph);
534 new (&hostJacTpetraBlocks(row,col)) LocalMatrixType(unmanagedMatrix);
537 hostBlockExistsInJac(row,col) = 1;
540 hostBlockExistsInJac(row,col) = 0;
544 typename PHX::View<LocalMatrixType**>
545 jacTpetraBlocks(
"panzer::ScatterResidual_BlockedTpetra<Jacobian>::jacTpetraBlocks",numFieldBlocks,numFieldBlocks);
546 Kokkos::deep_copy(jacTpetraBlocks,hostJacTpetraBlocks);
547 Kokkos::deep_copy(blockExistsInJac,hostBlockExistsInJac);
554 const auto& globalIndexers = globalIndexer_->getFieldDOFManagers();
555 auto blockOffsets_h = Kokkos::create_mirror_view(blockOffsets_);
556 Kokkos::deep_copy(blockOffsets_h, blockOffsets_);
557 for (
size_t block=0; block < globalIndexers.size(); ++block) {
558 const auto subviewOfBlockLIDs = Kokkos::subview(worksetLIDs_,Kokkos::ALL(), std::make_pair(blockOffsets_h(block),blockOffsets_h(block+1)));
559 globalIndexers[block]->getElementLIDs(localCellIds,subviewOfBlockLIDs);
563 for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
565 const int blockRowIndex = productVectorBlockIndex_[fieldIndex];
566 typename Tpetra::Vector<double,LO,GO,PHX::Device>::dual_view_type::t_dev kokkosResidual;
568 auto& tpetraResidual = *((rcp_dynamic_cast<Thyra::TpetraVector<RealType,LO,GO,NodeT>>(thyraBlockResidual->getNonconstVectorBlock(blockRowIndex),
true))->getTpetraVector());
569 kokkosResidual = tpetraResidual.getLocalViewDevice(Tpetra::Access::ReadWrite);
573 auto& tpetraDirichletCounter = *((rcp_dynamic_cast<Thyra::TpetraVector<RealType,LO,GO,NodeT>>(dirichletCounter_->getNonconstVectorBlock(blockRowIndex),
true))->getTpetraVector());
574 const auto& kokkosDirichletCounter = tpetraDirichletCounter.getLocalViewDevice(Tpetra::Access::ReadWrite);
577 const PHX::View<const int*> fieldOffsets = fieldOffsets_[fieldIndex];
578 const auto& basisIndices = basisIndexForMDFieldOffsets_[fieldIndex];
579 const PHX::View<const LO**> worksetLIDs = worksetLIDs_;
580 const PHX::View<const ScalarT**> fieldValues = scatterFields_[fieldIndex].get_static_view();
581 const PHX::View<const LO*> blockOffsets = blockOffsets_;
582 const auto& applyBC = applyBC_[fieldIndex].get_static_view();
583 const bool checkApplyBC = checkApplyBC_;
585 Kokkos::parallel_for(Kokkos::RangePolicy<PHX::Device>(0,workset.num_cells), KOKKOS_LAMBDA (
const int& cell) {
586 LO cLIDs[maxDerivativeArraySize_];
587 typename Sacado::ScalarType<ScalarT>::type
vals[maxDerivativeArraySize_];
589 for(
int basis=0; basis < static_cast<int>(fieldOffsets.size()); ++basis) {
590 const int block_row_offset = blockOffsets(blockRowIndex);
591 const int rowLID = worksetLIDs(cell,block_row_offset+fieldOffsets(basis));
596 const int basisIndex = basisIndices(basis);
600 if (!applyBC(cell,basisIndex))
604 typename FieldType::array_type::reference_type tmpFieldVal = fieldValues(cell,basisIndex);
607 kokkosResidual(rowLID,0) = tmpFieldVal.val();
609 kokkosDirichletCounter(rowLID,0) = 1.0;
612 for (
int blockColIndex=0; blockColIndex < numFieldBlocks; ++blockColIndex) {
613 if (blockExistsInJac(blockRowIndex,blockColIndex)) {
614 const auto& rowEntries = jacTpetraBlocks(blockRowIndex,blockColIndex).row(rowLID);
615 for (
int i=0; i < rowEntries.length; ++i)
621 for (
int blockColIndex=0; blockColIndex < numFieldBlocks; ++blockColIndex) {
622 if (blockExistsInJac(blockRowIndex,blockColIndex)) {
623 const int start = blockOffsets(blockColIndex);
624 const int stop = blockOffsets(blockColIndex+1);
625 const int sensSize = stop-start;
627 for (
int i=0; i < sensSize; ++i) {
628 cLIDs[i] = worksetLIDs(cell,start+i);
629 vals[i] = tmpFieldVal.fastAccessDx(start+i);
631 jacTpetraBlocks(blockRowIndex,blockColIndex).replaceValues(rowLID,cLIDs,sensSize,vals,
true,
true);
640 for (
int row=0; row < numFieldBlocks; ++row) {
641 for (
int col=0; col < numFieldBlocks; ++col) {
642 if (hostBlockExistsInJac(row,col)) {
643 hostJacTpetraBlocks(row,col).~CrsMatrix();
T & get(const std::string &name, T def_value)
bool nonnull(const std::shared_ptr< T > &p)
bool is_null(const std::shared_ptr< T > &p)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
ScatterDirichletResidual_BlockedTpetra(const Teuchos::RCP< const BlockedDOFManager > &)
FieldType
The type of discretization to use for a field pattern.
PHX::View< const int * > offsets
bool isParameter(const std::string &name) const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Kokkos::View< typename Sacado::ScalarType< ScalarT >::type **, Kokkos::LayoutRight, PHX::Device > vals
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()
bool isType(const std::string &name) const
#define TEUCHOS_ASSERT(assertion_test)
void postRegistrationSetup(typename TRAITS::SetupData, PHX::FieldManager< TRAITS > &)