43 #ifndef PANZER_SCATTER_RESIDUAL_TPETRA_IMPL_HPP 44 #define PANZER_SCATTER_RESIDUAL_TPETRA_IMPL_HPP 47 #include "Teuchos_Assert.hpp" 49 #include "Phalanx_DataLayout.hpp" 57 #include "Phalanx_DataLayout_MDALayout.hpp" 59 #include "Teuchos_FancyOStream.hpp" 61 #include "Tpetra_Vector.hpp" 62 #include "Tpetra_Map.hpp" 63 #include "Tpetra_CrsMatrix.hpp" 69 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
73 : globalIndexer_(indexer)
74 , globalDataKey_(
"Residual Scatter Container")
76 std::string scatterName = p.
get<std::string>(
"Scatter Name");
81 const std::vector<std::string>& names =
92 for (std::size_t eq = 0; eq < names.size(); ++eq) {
93 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
100 this->addEvaluatedField(*scatterHolder_);
102 if (p.
isType<std::string>(
"Global Data Key"))
103 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
105 this->setName(scatterName+
" Scatter Residual");
109 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
118 std::string fieldName = fieldMap_->find(
scatterFields_[fd].fieldTag().name())->second;
119 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
127 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
134 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(d.gedc.getDataObject(globalDataKey_));
136 if(tpetraContainer_==Teuchos::null) {
139 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(loc);
144 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
150 std::vector<LO> LIDs;
153 std::string blockId = this->wda(workset).block_id;
154 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
165 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
166 std::size_t cellLocalId = localCellIds[worksetCellIndex];
168 LIDs = globalIndexer_->getElementLIDs(cellLocalId);
171 for (std::size_t fieldIndex = 0; fieldIndex <
scatterFields_.size(); fieldIndex++) {
172 int fieldNum = fieldIds_[fieldIndex];
173 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
189 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
193 : globalIndexer_(indexer)
194 , globalDataKey_(
"Residual Scatter Container")
196 std::string scatterName = p.
get<std::string>(
"Scatter Name");
201 const std::vector<std::string>& names =
212 for (std::size_t eq = 0; eq < names.size(); ++eq) {
213 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
220 this->addEvaluatedField(*scatterHolder_);
222 if (p.
isType<std::string>(
"Global Data Key"))
223 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
225 this->setName(scatterName+
" Scatter Tangent");
229 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
238 std::string fieldName = fieldMap_->find(
scatterFields_[fd].fieldTag().name())->second;
239 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
247 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
252 using Teuchos::rcp_dynamic_cast;
257 std::vector<std::string> activeParameters =
260 dfdp_vectors_.clear();
261 for(std::size_t i=0;i<activeParameters.size();i++) {
263 rcp_dynamic_cast<LOC>(d.gedc.getDataObject(activeParameters[i]),
true)->get_f();
265 dfdp_vectors_.push_back(vec_array);
270 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
274 std::vector<LO> LIDs;
277 std::string blockId = this->wda(workset).block_id;
278 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
286 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
287 std::size_t cellLocalId = localCellIds[worksetCellIndex];
289 LIDs = globalIndexer_->getElementLIDs(cellLocalId);
292 for (std::size_t fieldIndex = 0; fieldIndex <
scatterFields_.size(); fieldIndex++) {
293 int fieldNum = fieldIds_[fieldIndex];
294 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
301 for(
int d=0;d<
value.size();d++)
302 dfdp_vectors_[d][lid] +=
value.fastAccessDx(d);
312 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
316 : globalIndexer_(indexer)
317 , globalDataKey_(
"Residual Scatter Container")
319 std::string scatterName = p.
get<std::string>(
"Scatter Name");
324 const std::vector<std::string>& names =
335 scratch_offsets_.resize(names.size());
336 for (std::size_t eq = 0; eq < names.size(); ++eq) {
337 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
344 this->addEvaluatedField(*scatterHolder_);
346 if (p.
isType<std::string>(
"Global Data Key"))
347 globalDataKey_ = p.
get<std::string>(
"Global Data Key");
349 this->setName(scatterName+
" Scatter Residual (Jacobian)");
353 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
360 const Workset & workset_0 = (*d.worksets_)[0];
361 std::string blockId = this->wda(workset_0).
block_id;
366 std::string fieldName = fieldMap_->find(
scatterFields_[fd].fieldTag().name())->second;
367 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
372 int fieldNum = fieldIds_[fd];
373 const std::vector<int> &
offsets = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
374 scratch_offsets_[fd] = Kokkos::View<int*,PHX::Device>(
"offsets",
offsets.size());
375 for(std::size_t i=0;i<
offsets.size();i++)
376 scratch_offsets_[fd](i) =
offsets[i];
379 scratch_lids_ = Kokkos::View<LO**,PHX::Device>(
"lids",
scatterFields_[0].dimension_0(),
380 globalIndexer_->getElementBlockGIDCount(blockId));
384 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
391 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(d.gedc.getDataObject(globalDataKey_));
393 if(tpetraContainer_==Teuchos::null) {
396 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(loc);
405 template <
typename ScalarT,
typename LO,
typename GO,
typename NodeT,
typename LocalMatrixT>
406 class ScatterResidual_Jacobian_Functor {
408 typedef typename PHX::Device execution_space;
409 typedef PHX::MDField<const ScalarT,Cell,NODE> FieldType;
412 Kokkos::View<double**, Kokkos::LayoutLeft,PHX::Device>
r_data;
415 Kokkos::View<const LO**,PHX::Device>
lids;
420 KOKKOS_INLINE_FUNCTION
421 void operator()(
const unsigned int cell)
const 424 typename Sacado::ScalarType<ScalarT>::type
vals[256];
425 int numIds =
lids.dimension_1();
427 for(
int i=0;i<numIds;i++)
428 cLIDs[i] =
lids(cell,i);
432 typename FieldType::array_type::reference_type scatterField =
field(cell,
basis);
438 Kokkos::atomic_add(&
r_data(lid,0), scatterField.val());
441 for(
int sensIndex=0;sensIndex<numIds;++sensIndex)
442 vals[sensIndex] = scatterField.fastAccessDx(sensIndex);
445 jac.sumIntoValues(lid, cLIDs,numIds,
vals,
true,
true);
453 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
460 std::vector<GO> GIDs;
461 std::vector<LO> cLIDs, rLIDs;
462 std::vector<double> jacRow;
465 std::string blockId = this->wda(workset).block_id;
466 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
477 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
478 std::size_t cellLocalId = localCellIds[worksetCellIndex];
481 rLIDs = globalIndexer_->getElementLIDs(cellLocalId);
485 for(std::size_t fieldIndex = 0; fieldIndex <
scatterFields_.size(); fieldIndex++) {
486 int fieldNum = fieldIds_[fieldIndex];
487 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
490 for(std::size_t rowBasisNum = 0; rowBasisNum < elmtOffset.size(); rowBasisNum++) {
492 int rowOffset = elmtOffset[rowBasisNum];
493 LO row = rLIDs[rowOffset];
497 r->sumIntoLocalValue(row,scatterField.val());
500 jacRow.resize(scatterField.size());
502 for(
int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
503 jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex);
506 Jac->sumIntoLocalValues(row, cLIDs, jacRow);
512 typedef typename LOC::CrsMatrixType::local_matrix_type LocalMatrixT;
515 std::string blockId = this->wda(workset).block_id;
520 globalIndexer_->getElementLIDs(this->wda(workset).cell_local_ids_k,scratch_lids_);
522 ScatterResidual_Jacobian_Functor<ScalarT,LO,GO,NodeT,LocalMatrixT> functor;
523 functor.fillResidual = (r!=Teuchos::null);
524 if(functor.fillResidual)
525 functor.r_data = r->template getLocalView<PHX::Device>();
526 functor.jac = Jac->getLocalMatrix();
527 functor.lids = scratch_lids_;
530 for(std::size_t fieldIndex = 0; fieldIndex <
scatterFields_.size(); fieldIndex++) {
531 functor.offsets = scratch_offsets_[fieldIndex];
534 Kokkos::parallel_for(workset.num_cells,functor);
panzer::Traits::Tangent::ScalarT ScalarT
std::vector< PHX::MDField< const ScalarT, panzer::Cell, panzer::Point > > scatterFields_
T & get(const std::string &name, T def_value)
void operator()(const FieldMultTag< NUM_FIELD_MULT > &, const std::size_t &cell) const
Kokkos::View< const LO **, PHX::Device > lids
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
bool isType(const std::string &name) const
Kokkos::View< double **, Kokkos::LayoutLeft, PHX::Device > r_data
Pushes residual values into the residual vector for a Newton-based solve.
Teuchos::RCP< const panzer::PureBasis > basis
Interpolates basis DOF values to IP DOF values.
panzer::Traits::Jacobian::ScalarT ScalarT
Kokkos::View< const int *, PHX::Device > offsets