43 #ifndef PANZER_SCATTER_RESIDUAL_TPETRA_IMPL_HPP 44 #define PANZER_SCATTER_RESIDUAL_TPETRA_IMPL_HPP 46 #include "Teuchos_RCP.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>
72 const Teuchos::ParameterList& p)
73 : globalIndexer_(indexer)
74 , globalDataKey_(
"Residual Scatter Container")
76 std::string scatterName = p.get<std::string>(
"Scatter Name");
78 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
81 const std::vector<std::string>& names =
82 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
85 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
87 Teuchos::RCP<PHX::DataLayout> dl =
88 p.get< Teuchos::RCP<const panzer::PureBasis> >(
"Basis")->functional;
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) {
138 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<
LOCPair_GlobalEvaluationData>(d.gedc.getDataObject(globalDataKey_),
true)->getGhostedLOC();
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;
156 Teuchos::RCP<typename LOC::VectorType> r = tpetraContainer_->
get_f();
157 Teuchos::ArrayRCP<double> r_array = r->get1dViewNonConst();
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>
192 const Teuchos::ParameterList& p)
193 : globalIndexer_(indexer)
194 , globalDataKey_(
"Residual Scatter Container")
196 std::string scatterName = p.get<std::string>(
"Scatter Name");
198 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
201 const std::vector<std::string>& names =
202 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
205 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
207 Teuchos::RCP<PHX::DataLayout> dl =
208 p.get< Teuchos::RCP<const panzer::PureBasis> >(
"Basis")->functional;
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++) {
262 RCP<typename LOC::VectorType> vec =
263 rcp_dynamic_cast<LOC>(d.gedc.getDataObject(activeParameters[i]),
true)->get_f();
264 Teuchos::ArrayRCP<double> vec_array = vec->get1dViewNonConst();
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>
315 const Teuchos::ParameterList& p)
316 : globalIndexer_(indexer)
317 , globalDataKey_(
"Residual Scatter Container")
319 std::string scatterName = p.get<std::string>(
"Scatter Name");
321 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
324 const std::vector<std::string>& names =
325 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
328 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
330 Teuchos::RCP<PHX::DataLayout> dl =
331 p.get< Teuchos::RCP<const panzer::PureBasis> >(
"Basis")->functional;
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) {
395 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<
LOCPair_GlobalEvaluationData>(d.gedc.getDataObject(globalDataKey_),
true)->getGhostedLOC();
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;
468 Teuchos::RCP<typename LOC::VectorType> r = tpetraContainer_->
get_f();
469 Teuchos::RCP<typename LOC::CrsMatrixType> Jac = tpetraContainer_->get_A();
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;
517 Teuchos::RCP<typename LOC::VectorType> r = tpetraContainer_->get_f();
518 Teuchos::RCP<typename LOC::CrsMatrixType> Jac = tpetraContainer_->get_A();
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_
void operator()(const FieldMultTag< NUM_FIELD_MULT > &, const std::size_t &cell) const
Kokkos::View< const LO **, PHX::Device > lids
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
const Teuchos::RCP< VectorType > get_f() const
Kokkos::View< const int *, PHX::Device > offsets