43 #ifndef PANZER_SCATTER_DIRICHLET_RESIDUAL_BLOCEDEPETRA_IMPL_HPP
44 #define PANZER_SCATTER_DIRICHLET_RESIDUAL_BLOCEDEPETRA_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"
61 #include "Phalanx_DataLayout_MDALayout.hpp"
63 #include "Thyra_SpmdVectorBase.hpp"
64 #include "Thyra_ProductVectorBase.hpp"
65 #include "Thyra_BlockedLinearOpBase.hpp"
66 #include "Thyra_get_Epetra_Operator.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) {
90 PHX::MDField<const ScalarT,Cell,NODE> scatterField = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
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) {
139 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
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) {
149 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string(
"APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
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] = Kokkos::View<int*,PHX::Device>(
"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] = Kokkos::View<int*,PHX::Device>(
"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] = Kokkos::View<int*,PHX::Device>(
"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);
218 PHX::Device::fence();
224 worksetLIDs_ = Kokkos::View<LO**,PHX::Device>(
"ScatterResidual_BlockedTpetra(Residual):worksetLIDs",
225 scatterFields_[0].extent(0),
226 maxElementBlockGIDCount);
230 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
236 = Teuchos::rcp_dynamic_cast<
ContainerType>(d.gedc->getDataObject(
"Dirichlet Counter"),
true);
242 blockedContainer_ = Teuchos::rcp_dynamic_cast<
const ContainerType>(d.gedc->getDataObject(globalDataKey_),
true);
247 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
252 using Teuchos::rcp_dynamic_cast;
256 const auto& localCellIds = this->wda(workset).cell_local_ids_k;
259 rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer_->get_f(),
true) :
260 rcp_dynamic_cast<ProductVectorBase<double> >(blockedContainer_->get_x(),
true);
263 int currentWorksetLIDSubBlock = -1;
264 for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
266 if (productVectorBlockIndex_[fieldIndex] != currentWorksetLIDSubBlock) {
267 fieldGlobalIndexers_[fieldIndex]->getElementLIDs(localCellIds,worksetLIDs_);
268 currentWorksetLIDSubBlock = productVectorBlockIndex_[fieldIndex];
272 const auto& tpetraScatterTarget = *((rcp_dynamic_cast<Thyra::TpetraVector<RealType,LO,GO,NodeT>>(thyraScatterTarget->getNonconstVectorBlock(productVectorBlockIndex_[fieldIndex]),
true))->getTpetraVector());
273 const auto& kokkosScatterTarget = tpetraScatterTarget.template getLocalView<PHX::mem_space>();
276 const auto& tpetraDirichletCounter = *((rcp_dynamic_cast<Thyra::TpetraVector<RealType,LO,GO,NodeT>>(dirichletCounter_->getNonconstVectorBlock(productVectorBlockIndex_[fieldIndex]),
true))->getTpetraVector());
277 const auto& kokkosDirichletCounter = tpetraDirichletCounter.template getLocalView<PHX::mem_space>();
280 const auto& fieldOffsets = fieldOffsets_[fieldIndex];
281 const auto& basisIndices = basisIndexForMDFieldOffsets_[fieldIndex];
282 const auto& worksetLIDs = worksetLIDs_;
283 const auto& fieldValues = scatterFields_[fieldIndex];
284 const auto& applyBC = applyBC_[fieldIndex].get_static_view();
285 const bool checkApplyBC = checkApplyBC_;
289 Kokkos::parallel_for(Kokkos::RangePolicy<PHX::Device>(0,workset.num_cells), KOKKOS_LAMBDA (
const int& cell) {
290 for (
int basis=0; basis < static_cast<int>(fieldOffsets.size()); ++basis) {
291 const int lid = worksetLIDs(cell,fieldOffsets(basis));
294 const int basisIndex = basisIndices(basis);
298 if (!applyBC(cell,basisIndex))
301 kokkosScatterTarget(lid,0) = fieldValues(cell,basisIndex);
302 kokkosDirichletCounter(lid,0) = 1.0;
308 Kokkos::parallel_for(Kokkos::RangePolicy<PHX::Device>(0,workset.num_cells), KOKKOS_LAMBDA (
const int& cell) {
309 for (
int basis=0; basis < static_cast<int>(fieldOffsets.size()); ++basis) {
310 const int lid = worksetLIDs(cell,fieldOffsets(basis));
313 kokkosScatterTarget(lid,0) = fieldValues(cell,basis);
314 kokkosDirichletCounter(lid,0) = 1.0;
327 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
331 : globalIndexer_(indexer)
332 , globalDataKey_(
"Residual Scatter Container")
334 std::string scatterName = p.
get<std::string>(
"Scatter Name");
339 const std::vector<std::string>& names =
348 side_subcell_dim_ = p.
get<
int>(
"Side Subcell Dimension");
349 local_side_id_ = p.
get<
int>(
"Local Side ID");
352 scatterFields_.resize(names.size());
353 for (std::size_t eq = 0; eq < names.size(); ++eq) {
354 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
357 this->addDependentField(scatterFields_[eq]);
360 checkApplyBC_ = p.
get<
bool>(
"Check Apply BC");
361 applyBC_.resize(names.size());
363 for (std::size_t eq = 0; eq < names.size(); ++eq) {
364 applyBC_[eq] = PHX::MDField<const bool,Cell,NODE>(std::string(
"APPLY_BC_")+fieldMap_->find(names[eq])->second,dl);
365 this->addDependentField(applyBC_[eq]);
370 this->addEvaluatedField(*scatterHolder_);
372 if (p.
isType<std::string>(
"Global Data Key"))
373 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
375 this->setName(scatterName+
" Scatter Residual (Jacobian)");
379 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
384 const Workset & workset_0 = (*d.worksets_)[0];
385 const std::string blockId = this->wda(workset_0).
block_id;
387 fieldIds_.resize(scatterFields_.size());
388 fieldOffsets_.resize(scatterFields_.size());
389 basisIndexForMDFieldOffsets_.resize(scatterFields_.size());
390 productVectorBlockIndex_.resize(scatterFields_.size());
391 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
393 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
394 const int globalFieldNum = globalIndexer_->getFieldNum(fieldName);
395 productVectorBlockIndex_[fd] = globalIndexer_->getFieldBlock(globalFieldNum);
396 const auto& fieldGlobalIndexer = globalIndexer_->getFieldDOFManagers()[productVectorBlockIndex_[fd]];
397 fieldIds_[fd] = fieldGlobalIndexer->getFieldNum(fieldName);
399 const auto& offsetPair = fieldGlobalIndexer->getGIDFieldOffsets_closure(blockId,fieldIds_[fd],side_subcell_dim_,local_side_id_);
401 const auto&
offsets = offsetPair.first;
402 fieldOffsets_[fd] = Kokkos::View<int*,PHX::Device>(
"ScatterDirichletResidual_BlockedTpetra(Jacobian):fieldOffsets",
offsets.size());
403 auto hostOffsets = Kokkos::create_mirror_view(fieldOffsets_[fd]);
404 for (std::size_t i=0; i <
offsets.size(); ++i)
406 Kokkos::deep_copy(fieldOffsets_[fd], hostOffsets);
409 const auto& basisIndex = offsetPair.second;
410 basisIndexForMDFieldOffsets_[fd] = Kokkos::View<int*,PHX::Device>(
"ScatterDirichletResidual_BlockedTpetra(Jacobian):basisIndexForMDFieldOffsets",basisIndex.size());
411 auto hostBasisIndex = Kokkos::create_mirror_view(basisIndexForMDFieldOffsets_[fd]);
412 for (std::size_t i=0; i < basisIndex.size(); ++i)
413 hostBasisIndex(i) = basisIndex[i];
414 Kokkos::deep_copy(basisIndexForMDFieldOffsets_[fd], hostBasisIndex);
421 int elementBlockGIDCount = 0;
422 for (
const auto blockDOFMgr : globalIndexer_->getFieldDOFManagers())
423 elementBlockGIDCount += blockDOFMgr->getElementBlockGIDCount(blockId);
425 worksetLIDs_ = Kokkos::View<LO**,PHX::Device>(
"ScatterDirichletResidual_BlockedTpetra(Jacobian):worksetLIDs",
426 scatterFields_[0].extent(0),
427 elementBlockGIDCount);
430 const auto& blockGlobalIndexers = globalIndexer_->getFieldDOFManagers();
431 const int numBlocks =
static_cast<int>(globalIndexer_->getFieldDOFManagers().size());
432 blockOffsets_ = Kokkos::View<LO*,PHX::Device>(
"ScatterDirichletResidual_BlockedTpetra(Jacobian):blockOffsets_",
434 const auto hostBlockOffsets = Kokkos::create_mirror_view(blockOffsets_);
435 for (
int blk=0;blk<numBlocks;blk++) {
436 int blockOffset = globalIndexer_->getBlockGIDOffset(blockId,blk);
437 hostBlockOffsets(blk) = blockOffset;
439 blockOffsets_(numBlocks) = blockOffsets_(numBlocks-1) + blockGlobalIndexers[blockGlobalIndexers.size()-1]->getElementBlockGIDCount(blockId);
440 Kokkos::deep_copy(blockOffsets_,hostBlockOffsets);
445 for (
int blk=0;blk<numBlocks;blk++) {
446 const int blockDerivativeSize = hostBlockOffsets(blk+1) - hostBlockOffsets(blk);
448 "ERROR: the derivative dimension for sub block "
449 << blk <<
"with a value of " << blockDerivativeSize
450 <<
"is larger than the size allocated for cLIDs and vals "
451 <<
"in the evaluate call! You must manually increase the "
452 <<
"size and recompile!");
455 PHX::Device::fence();
459 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
465 = Teuchos::rcp_dynamic_cast<
const ContainerType>(d.gedc->getDataObject(
"Dirichlet Counter"),
true);
471 blockedContainer_ = Teuchos::rcp_dynamic_cast<
const ContainerType>(d.gedc->getDataObject(globalDataKey_),
true);
476 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
481 using Teuchos::rcp_dynamic_cast;
486 const auto& localCellIds = this->wda(workset).cell_local_ids_k;
488 const int numFieldBlocks = globalIndexer_->getNumFieldBlocks();
498 using LocalMatrixType = KokkosSparse::CrsMatrix<double,LO,PHX::Device,Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
499 typename Kokkos::View<LocalMatrixType**,PHX::Device>::HostMirror
500 hostJacTpetraBlocks(
"panzer::ScatterResidual_BlockTpetra<Jacobian>::hostJacTpetraBlocks", numFieldBlocks,numFieldBlocks);
502 Kokkos::View<int**,PHX::Device> blockExistsInJac = Kokkos::View<int**,PHX::Device>(
"blockExistsInJac_",numFieldBlocks,numFieldBlocks);
503 auto hostBlockExistsInJac = Kokkos::create_mirror_view(blockExistsInJac);
505 for (
int row=0; row < numFieldBlocks; ++row) {
506 for (
int col=0; col < numFieldBlocks; ++col) {
507 const auto thyraTpetraOperator = rcp_dynamic_cast<Thyra::TpetraLinearOp<double,LO,GO,NodeT>>(Jac->getNonconstBlock(row,col),
false);
508 if (
nonnull(thyraTpetraOperator)) {
524 const auto tpetraCrsMatrix = rcp_dynamic_cast<Tpetra::CrsMatrix<double,LO,GO,NodeT>>(thyraTpetraOperator->getTpetraOperator(),
true);
525 const auto managedMatrix = tpetraCrsMatrix->getLocalMatrix();
526 const auto managedGraph = managedMatrix.graph;
529 using StaticCrsGraphType =
typename LocalMatrixType::StaticCrsGraphType;
530 StaticCrsGraphType unmanagedGraph;
531 unmanagedGraph.entries =
typename StaticCrsGraphType::entries_type(managedGraph.entries.data(),managedGraph.entries.extent(0));
532 unmanagedGraph.row_map =
typename StaticCrsGraphType::row_map_type(managedGraph.row_map.data(),managedGraph.row_map.extent(0));
533 unmanagedGraph.row_block_offsets =
typename StaticCrsGraphType::row_block_type(managedGraph.row_block_offsets.data(),managedGraph.row_block_offsets.extent(0));
535 typename LocalMatrixType::values_type unmanagedValues(managedMatrix.values.data(),managedMatrix.values.extent(0));
536 LocalMatrixType unmanagedMatrix(managedMatrix.values.label(), managedMatrix.numCols(), unmanagedValues, unmanagedGraph);
537 new (&hostJacTpetraBlocks(row,col)) LocalMatrixType(unmanagedMatrix);
540 hostBlockExistsInJac(row,col) = 1;
543 hostBlockExistsInJac(row,col) = 0;
547 typename Kokkos::View<LocalMatrixType**,PHX::Device>
548 jacTpetraBlocks(
"panzer::ScatterResidual_BlockedTpetra<Jacobian>::jacTpetraBlocks",numFieldBlocks,numFieldBlocks);
549 Kokkos::deep_copy(jacTpetraBlocks,hostJacTpetraBlocks);
550 Kokkos::deep_copy(blockExistsInJac,hostBlockExistsInJac);
551 PHX::Device::fence();
558 const auto& globalIndexers = globalIndexer_->getFieldDOFManagers();
559 for (
size_t block=0; block < globalIndexers.size(); ++block) {
560 const auto subviewOfBlockLIDs = Kokkos::subview(worksetLIDs_,Kokkos::ALL(), std::make_pair(blockOffsets_(block),blockOffsets_(block+1)));
561 globalIndexers[block]->getElementLIDs(localCellIds,subviewOfBlockLIDs);
565 for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
567 const int blockRowIndex = productVectorBlockIndex_[fieldIndex];
568 typename Tpetra::Vector<double,LO,GO,PHX::Device>::dual_view_type::t_dev kokkosResidual;
570 const auto& tpetraResidual = *((rcp_dynamic_cast<Thyra::TpetraVector<RealType,LO,GO,NodeT>>(thyraBlockResidual->getNonconstVectorBlock(blockRowIndex),
true))->getTpetraVector());
571 kokkosResidual = tpetraResidual.template getLocalView<PHX::mem_space>();
575 const auto& tpetraDirichletCounter = *((rcp_dynamic_cast<Thyra::TpetraVector<RealType,LO,GO,NodeT>>(dirichletCounter_->getNonconstVectorBlock(blockRowIndex),
true))->getTpetraVector());
576 const auto& kokkosDirichletCounter = tpetraDirichletCounter.template getLocalView<PHX::mem_space>();
579 const Kokkos::View<const int*,PHX::Device> fieldOffsets = fieldOffsets_[fieldIndex];
580 const auto& basisIndices = basisIndexForMDFieldOffsets_[fieldIndex];
581 const Kokkos::View<const LO**,PHX::Device> worksetLIDs = worksetLIDs_;
582 const PHX::View<const ScalarT**> fieldValues = scatterFields_[fieldIndex].get_static_view();
583 const Kokkos::View<const LO*,PHX::Device> blockOffsets = blockOffsets_;
584 const auto& applyBC = applyBC_[fieldIndex].get_static_view();
585 const bool checkApplyBC = checkApplyBC_;
587 Kokkos::parallel_for(Kokkos::RangePolicy<PHX::Device>(0,workset.num_cells), KOKKOS_LAMBDA (
const int& cell) {
588 LO cLIDs[maxDerivativeArraySize_];
589 typename Sacado::ScalarType<ScalarT>::type vals[maxDerivativeArraySize_];
591 for(
int basis=0; basis < static_cast<int>(fieldOffsets.size()); ++basis) {
592 const int rowLID = worksetLIDs(cell,fieldOffsets(basis));
597 const int basisIndex = basisIndices(basis);
601 if (!applyBC(cell,basisIndex))
604 typedef PHX::MDField<const ScalarT,Cell,NODE>
FieldType;
605 typename FieldType::array_type::reference_type tmpFieldVal = fieldValues(cell,basisIndex);
608 kokkosResidual(rowLID,0) = tmpFieldVal.val();
610 kokkosDirichletCounter(rowLID,0) = 1.0;
613 for (
int blockColIndex=0; blockColIndex < numFieldBlocks; ++blockColIndex) {
614 if (blockExistsInJac(blockRowIndex,blockColIndex)) {
615 const auto& rowEntries = jacTpetraBlocks(blockRowIndex,blockColIndex).row(rowLID);
616 for (
int i=0; i < rowEntries.length; ++i)
622 for (
int blockColIndex=0; blockColIndex < numFieldBlocks; ++blockColIndex) {
623 if (blockExistsInJac(blockRowIndex,blockColIndex)) {
624 const int start = blockOffsets(blockColIndex);
625 const int stop = blockOffsets(blockColIndex+1);
626 const int sensSize = stop-start;
628 for (
int i=0; i < sensSize; ++i) {
629 cLIDs[i] = worksetLIDs(cell,start+i);
630 vals[i] = tmpFieldVal.fastAccessDx(start+i);
632 jacTpetraBlocks(blockRowIndex,blockColIndex).replaceValues(rowLID,cLIDs,sensSize,vals,
true,
true);
641 for (
int row=0; row < numFieldBlocks; ++row) {
642 for (
int col=0; col < numFieldBlocks; ++col) {
643 if (hostBlockExistsInJac(row,col)) {
644 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)
FieldType
The type of discretization to use for a field pattern.
bool isParameter(const std::string &name) const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void evaluateFields(typename TRAITS::EvalData)
bool nonnull(const boost::shared_ptr< T > &p)
bool isType(const std::string &name) const
Pushes residual values into the residual vector for a Newton-based solve.
#define TEUCHOS_ASSERT(assertion_test)
void postRegistrationSetup(typename TRAITS::SetupData, PHX::FieldManager< TRAITS > &)
Kokkos::View< const int *, PHX::Device > offsets
ScatterDirichletResidual_BlockedTpetra(const Teuchos::RCP< const BlockedDOFManager< LO, GO > > &)