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) {
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>
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) {
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] = Kokkos::View<int*,PHX::Device>(
"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);
172 PHX::Device::fence();
178 worksetLIDs_ = Kokkos::View<LO**,PHX::Device>(
"ScatterResidual_BlockedTpetra(Residual):worksetLIDs",
179 scatterFields_[0].extent(0),
180 maxElementBlockGIDCount);
184 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
189 blockedContainer_ = Teuchos::rcp_dynamic_cast<
const ContainerType>(d.gedc->getDataObject(globalDataKey_),
true);
193 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
198 using Teuchos::rcp_dynamic_cast;
202 const auto& localCellIds = this->wda(workset).cell_local_ids_k;
206 int currentWorksetLIDSubBlock = -1;
207 for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
209 if (productVectorBlockIndex_[fieldIndex] != currentWorksetLIDSubBlock) {
210 fieldGlobalIndexers_[fieldIndex]->getElementLIDs(localCellIds,worksetLIDs_);
211 currentWorksetLIDSubBlock = productVectorBlockIndex_[fieldIndex];
214 const auto& tpetraResidual = *((rcp_dynamic_cast<Thyra::TpetraVector<RealType,LO,GO,NodeT>>(thyraBlockResidual->getNonconstVectorBlock(productVectorBlockIndex_[fieldIndex]),
true))->getTpetraVector());
215 const auto& kokkosResidual = tpetraResidual.template getLocalView<PHX::mem_space>();
218 const auto& fieldOffsets = fieldOffsets_[fieldIndex];
219 const auto& worksetLIDs = worksetLIDs_;
220 const auto& fieldValues = scatterFields_[fieldIndex];
222 Kokkos::parallel_for(Kokkos::RangePolicy<PHX::Device>(0,workset.num_cells), KOKKOS_LAMBDA (
const int& cell) {
223 for(
int basis=0; basis < static_cast<int>(fieldOffsets.size()); ++basis) {
224 const int lid = worksetLIDs(cell,fieldOffsets(basis));
225 Kokkos::atomic_add(&kokkosResidual(lid,0), fieldValues(cell,basis));
235 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
239 : globalIndexer_(indexer)
240 , globalDataKey_(
"Residual Scatter Container")
242 std::string scatterName = p.
get<std::string>(
"Scatter Name");
247 const std::vector<std::string>& names =
257 scatterFields_.resize(names.size());
258 for (std::size_t eq = 0; eq < names.size(); ++eq) {
259 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
262 this->addDependentField(scatterFields_[eq]);
266 this->addEvaluatedField(*scatterHolder_);
268 if (p.
isType<std::string>(
"Global Data Key"))
269 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
271 this->setName(scatterName+
" Scatter Residual (Jacobian)");
275 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
280 const Workset & workset_0 = (*d.worksets_)[0];
281 const std::string blockId = this->wda(workset_0).
block_id;
283 fieldIds_.resize(scatterFields_.size());
284 fieldOffsets_.resize(scatterFields_.size());
285 productVectorBlockIndex_.resize(scatterFields_.size());
286 for (std::size_t fd=0; fd < scatterFields_.size(); ++fd) {
287 const std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
288 const int globalFieldNum = globalIndexer_->getFieldNum(fieldName);
289 productVectorBlockIndex_[fd] = globalIndexer_->getFieldBlock(globalFieldNum);
290 const auto& fieldGlobalIndexer = globalIndexer_->getFieldDOFManagers()[productVectorBlockIndex_[fd]];
291 fieldIds_[fd] = fieldGlobalIndexer->getFieldNum(fieldName);
293 const std::vector<int>&
offsets = globalIndexer_->getGIDFieldOffsets(blockId,globalFieldNum);
294 fieldOffsets_[fd] = Kokkos::View<int*,PHX::Device>(
"ScatterResidual_BlockedTpetra(Jacobian):fieldOffsets",offsets.size());
295 auto hostOffsets = Kokkos::create_mirror_view(fieldOffsets_[fd]);
296 for (std::size_t i=0; i < offsets.size(); ++i)
297 hostOffsets(i) = offsets[i];
298 Kokkos::deep_copy(fieldOffsets_[fd], hostOffsets);
299 PHX::Device::fence();
305 int elementBlockGIDCount = 0;
306 for (
const auto blockDOFMgr : globalIndexer_->getFieldDOFManagers())
307 elementBlockGIDCount += blockDOFMgr->getElementBlockGIDCount(blockId);
309 worksetLIDs_ = Kokkos::View<LO**,PHX::Device>(
"ScatterResidual_BlockedTpetra(Jacobian):worksetLIDs",
310 scatterFields_[0].extent(0),
311 elementBlockGIDCount);
314 const auto& blockGlobalIndexers = globalIndexer_->getFieldDOFManagers();
315 const int numBlocks =
static_cast<int>(globalIndexer_->getFieldDOFManagers().size());
316 blockOffsets_ = Kokkos::View<LO*,PHX::Device>(
"ScatterResidual_BlockedTpetra(Jacobian):blockOffsets_",
318 const auto hostBlockOffsets = Kokkos::create_mirror_view(blockOffsets_);
319 for (
int blk=0;blk<numBlocks;blk++) {
320 int blockOffset = globalIndexer_->getBlockGIDOffset(blockId,blk);
321 hostBlockOffsets(blk) = blockOffset;
323 blockOffsets_(numBlocks) = blockOffsets_(numBlocks-1) + blockGlobalIndexers[blockGlobalIndexers.size()-1]->getElementBlockGIDCount(blockId);
324 Kokkos::deep_copy(blockOffsets_,hostBlockOffsets);
329 for (
int blk=0;blk<numBlocks;blk++) {
330 const int blockDerivativeSize = hostBlockOffsets(blk+1) - hostBlockOffsets(blk);
332 "ERROR: the derivative dimension for sub block "
333 << blk <<
"with a value of " << blockDerivativeSize
334 <<
"is larger than the size allocated for cLIDs and vals "
335 <<
"in the evaluate call! You must manually increase the "
336 <<
"size and recompile!");
339 PHX::Device::fence();
343 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
348 using Teuchos::rcp_dynamic_cast;
351 blockedContainer_ = rcp_dynamic_cast<
const ContainerType>(d.gedc->getDataObject(globalDataKey_));
353 if(blockedContainer_==Teuchos::null) {
360 template <
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
365 using Teuchos::rcp_dynamic_cast;
370 const auto& localCellIds = this->wda(workset).cell_local_ids_k;
372 const int numFieldBlocks = globalIndexer_->getNumFieldBlocks();
382 using LocalMatrixType = KokkosSparse::CrsMatrix<double,LO,PHX::Device,Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
383 typename Kokkos::View<LocalMatrixType**,PHX::Device>::HostMirror
384 hostJacTpetraBlocks(
"panzer::ScatterResidual_BlockTpetra<Jacobian>::hostJacTpetraBlocks", numFieldBlocks,numFieldBlocks);
386 Kokkos::View<int**,PHX::Device> blockExistsInJac = Kokkos::View<int**,PHX::Device>(
"blockExistsInJac_",numFieldBlocks,numFieldBlocks);
387 auto hostBlockExistsInJac = Kokkos::create_mirror_view(blockExistsInJac);
389 for (
int row=0; row < numFieldBlocks; ++row) {
390 for (
int col=0; col < numFieldBlocks; ++col) {
391 const auto thyraTpetraOperator = rcp_dynamic_cast<Thyra::TpetraLinearOp<double,LO,GO,NodeT>>(Jac->getNonconstBlock(row,col),
false);
392 if (
nonnull(thyraTpetraOperator)) {
408 const auto tpetraCrsMatrix = rcp_dynamic_cast<Tpetra::CrsMatrix<double,LO,GO,NodeT>>(thyraTpetraOperator->getTpetraOperator(),
true);
409 const auto managedMatrix = tpetraCrsMatrix->getLocalMatrix();
410 const auto managedGraph = managedMatrix.graph;
413 using StaticCrsGraphType =
typename LocalMatrixType::StaticCrsGraphType;
414 StaticCrsGraphType unmanagedGraph;
415 unmanagedGraph.entries =
typename StaticCrsGraphType::entries_type(managedGraph.entries.data(),managedGraph.entries.extent(0));
416 unmanagedGraph.row_map =
typename StaticCrsGraphType::row_map_type(managedGraph.row_map.data(),managedGraph.row_map.extent(0));
417 unmanagedGraph.row_block_offsets =
typename StaticCrsGraphType::row_block_type(managedGraph.row_block_offsets.data(),managedGraph.row_block_offsets.extent(0));
419 typename LocalMatrixType::values_type unmanagedValues(managedMatrix.values.data(),managedMatrix.values.extent(0));
420 LocalMatrixType unmanagedMatrix(managedMatrix.values.label(), managedMatrix.numCols(), unmanagedValues, unmanagedGraph);
421 new (&hostJacTpetraBlocks(row,col)) LocalMatrixType(unmanagedMatrix);
424 hostBlockExistsInJac(row,col) = 1;
427 hostBlockExistsInJac(row,col) = 0;
431 typename Kokkos::View<LocalMatrixType**,PHX::Device>
432 jacTpetraBlocks(
"panzer::ScatterResidual_BlockedTpetra<Jacobian>::jacTpetraBlocks",numFieldBlocks,numFieldBlocks);
433 Kokkos::deep_copy(jacTpetraBlocks,hostJacTpetraBlocks);
434 Kokkos::deep_copy(blockExistsInJac,hostBlockExistsInJac);
435 PHX::Device::fence();
442 const auto& globalIndexers = globalIndexer_->getFieldDOFManagers();
443 for (
size_t block=0; block < globalIndexers.size(); ++block) {
444 const auto subviewOfBlockLIDs = Kokkos::subview(worksetLIDs_,Kokkos::ALL(), std::make_pair(blockOffsets_(block),blockOffsets_(block+1)));
445 globalIndexers[block]->getElementLIDs(localCellIds,subviewOfBlockLIDs);
449 for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
451 const int blockRowIndex = productVectorBlockIndex_[fieldIndex];
452 typename Tpetra::Vector<double,LO,GO,PHX::Device>::dual_view_type::t_dev kokkosResidual;
454 const auto& tpetraResidual = *((rcp_dynamic_cast<Thyra::TpetraVector<RealType,LO,GO,NodeT>>(thyraBlockResidual->getNonconstVectorBlock(blockRowIndex),
true))->getTpetraVector());
455 kokkosResidual = tpetraResidual.template getLocalView<PHX::mem_space>();
459 const Kokkos::View<const int*,PHX::Device> fieldOffsets = fieldOffsets_[fieldIndex];
460 const Kokkos::View<const LO**,PHX::Device> worksetLIDs = worksetLIDs_;
461 const PHX::View<const ScalarT**> fieldValues = scatterFields_[fieldIndex].get_static_view();
462 const Kokkos::View<const LO*,PHX::Device> blockOffsets = blockOffsets_;
464 Kokkos::parallel_for(Kokkos::RangePolicy<PHX::Device>(0,workset.num_cells), KOKKOS_LAMBDA (
const int& cell) {
466 typename Sacado::ScalarType<ScalarT>::type vals[256];
468 for(
int basis=0; basis < static_cast<int>(fieldOffsets.size()); ++basis) {
469 typedef PHX::MDField<const ScalarT,Cell,NODE>
FieldType;
470 typename FieldType::array_type::reference_type tmpFieldVal = fieldValues(cell,basis);
471 const int rowLID = worksetLIDs(cell,fieldOffsets(basis));
474 Kokkos::atomic_add(&kokkosResidual(rowLID,0), tmpFieldVal.val());
476 for (
int blockColIndex=0; blockColIndex < numFieldBlocks; ++blockColIndex) {
477 if (blockExistsInJac(blockRowIndex,blockColIndex)) {
478 const int start = blockOffsets(blockColIndex);
479 const int stop = blockOffsets(blockColIndex+1);
480 const int sensSize = stop-start;
482 for (
int i=0; i < sensSize; ++i) {
483 cLIDs[i] = worksetLIDs(cell,start+i);
484 vals[i] = tmpFieldVal.fastAccessDx(start+i);
486 jacTpetraBlocks(blockRowIndex,blockColIndex).sumIntoValues(rowLID,cLIDs,sensSize,vals,
true,
true);
495 for (
int row=0; row < numFieldBlocks; ++row) {
496 for (
int col=0; col < numFieldBlocks; ++col) {
497 if (hostBlockExistsInJac(row,col)) {
498 hostJacTpetraBlocks(row,col).~CrsMatrix();
Pushes residual values into the residual vector for a Newton-based solve.
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 > &)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
FieldType
The type of discretization to use for a field pattern.
Teuchos::RCP< LinearObjContainer > getGhostedLOC() 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)
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
ScatterResidual_BlockedTpetra(const Teuchos::RCP< const BlockedDOFManager< LO, GO > > &)
Kokkos::View< const int *, PHX::Device > offsets