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" 58 #include "Phalanx_DataLayout_MDALayout.hpp" 60 #include "Teuchos_FancyOStream.hpp" 62 #include "Tpetra_Vector.hpp" 63 #include "Tpetra_Map.hpp" 64 #include "Tpetra_CrsMatrix.hpp" 70 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
73 const Teuchos::ParameterList& p)
74 : globalIndexer_(indexer)
75 , globalDataKey_(
"Residual Scatter Container")
77 std::string scatterName = p.get<std::string>(
"Scatter Name");
79 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
82 const std::vector<std::string>& names =
83 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
86 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
88 Teuchos::RCP<PHX::DataLayout> dl =
89 p.get< Teuchos::RCP<const panzer::PureBasis> >(
"Basis")->functional;
92 scatterFields_.resize(names.size());
93 scratch_offsets_.resize(names.size());
94 for (std::size_t eq = 0; eq < names.size(); ++eq) {
95 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
98 this->addDependentField(scatterFields_[eq]);
102 this->addEvaluatedField(*scatterHolder_);
104 if (p.isType<std::string>(
"Global Data Key"))
105 globalDataKey_ = p.get<std::string>(
"Global Data Key");
107 this->setName(scatterName+
" Scatter Residual");
111 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
116 fieldIds_.resize(scatterFields_.size());
117 const Workset & workset_0 = (*d.worksets_)[0];
118 std::string blockId = this->wda(workset_0).
block_id;
122 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
124 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
125 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
127 const std::vector<int> &
offsets = globalIndexer_->getGIDFieldOffsets(blockId,fieldIds_[fd]);
128 scratch_offsets_[fd] = PHX::View<int*>(
"offsets",
offsets.size());
129 Kokkos::deep_copy(scratch_offsets_[fd], Kokkos::View<const int*, Kokkos::HostSpace, Kokkos::MemoryUnmanaged>(
offsets.data(),
offsets.size()));
131 scratch_lids_ = PHX::View<LO**>(
"lids",scatterFields_[0].extent(0),
132 globalIndexer_->getElementBlockGIDCount(blockId));
137 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
144 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(d.gedc->getDataObject(globalDataKey_));
146 if(tpetraContainer_==Teuchos::null) {
148 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<
LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),
true)->getGhostedLOC();
149 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(loc);
158 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
161 const Teuchos::ParameterList& p)
162 : globalIndexer_(indexer)
163 , globalDataKey_(
"Residual Scatter Container")
165 std::string scatterName = p.get<std::string>(
"Scatter Name");
167 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
170 const std::vector<std::string>& names =
171 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
174 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
176 Teuchos::RCP<PHX::DataLayout> dl =
177 p.get< Teuchos::RCP<const panzer::PureBasis> >(
"Basis")->functional;
180 scatterFields_.resize(names.size());
181 for (std::size_t eq = 0; eq < names.size(); ++eq) {
182 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
185 this->addDependentField(scatterFields_[eq]);
189 this->addEvaluatedField(*scatterHolder_);
191 if (p.isType<std::string>(
"Global Data Key"))
192 globalDataKey_ = p.get<std::string>(
"Global Data Key");
194 this->setName(scatterName+
" Scatter Tangent");
198 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
203 fieldIds_.resize(scatterFields_.size());
205 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
207 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
208 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
213 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
218 using Teuchos::rcp_dynamic_cast;
223 std::vector<std::string> activeParameters =
226 dfdp_vectors_.clear();
227 for(std::size_t i=0;i<activeParameters.size();i++) {
228 RCP<typename LOC::VectorType> vec =
229 rcp_dynamic_cast<LOC>(d.gedc->getDataObject(activeParameters[i]),
true)->get_f();
230 Teuchos::ArrayRCP<double> vec_array = vec->get1dViewNonConst();
231 dfdp_vectors_.push_back(vec_array);
236 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
241 std::string blockId = this->wda(workset).block_id;
242 const std::vector<std::size_t> & localCellIds = this->wda(workset).cell_local_ids;
250 for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
251 std::size_t cellLocalId = localCellIds[worksetCellIndex];
253 auto LIDs = globalIndexer_->getElementLIDs(cellLocalId);
256 for (std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
257 int fieldNum = fieldIds_[fieldIndex];
258 const std::vector<int> & elmtOffset = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
261 for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
262 int offset = elmtOffset[basis];
263 LO lid = LIDs[offset];
264 ScalarT value = (scatterFields_[fieldIndex])(worksetCellIndex,basis);
265 for(
int d=0;d<value.size();d++)
266 dfdp_vectors_[d][lid] += value.fastAccessDx(d);
276 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
279 const Teuchos::ParameterList& p)
280 : globalIndexer_(indexer)
281 , globalDataKey_(
"Residual Scatter Container")
282 , my_derivative_size_(0)
283 , other_derivative_size_(0)
285 std::string scatterName = p.get<std::string>(
"Scatter Name");
287 Teuchos::rcp(
new PHX::Tag<ScalarT>(scatterName,Teuchos::rcp(
new PHX::MDALayout<Dummy>(0))));
290 const std::vector<std::string>& names =
291 *(p.get< Teuchos::RCP< std::vector<std::string> > >(
"Dependent Names"));
294 fieldMap_ = p.get< Teuchos::RCP< std::map<std::string,std::string> > >(
"Dependent Map");
296 Teuchos::RCP<PHX::DataLayout> dl =
297 p.get< Teuchos::RCP<const panzer::PureBasis> >(
"Basis")->functional;
300 scatterFields_.resize(names.size());
301 scratch_offsets_.resize(names.size());
302 for (std::size_t eq = 0; eq < names.size(); ++eq) {
303 scatterFields_[eq] = PHX::MDField<const ScalarT,Cell,NODE>(names[eq],dl);
306 this->addDependentField(scatterFields_[eq]);
310 this->addEvaluatedField(*scatterHolder_);
312 if (p.isType<std::string>(
"Global Data Key"))
313 globalDataKey_ = p.get<std::string>(
"Global Data Key");
315 this->setName(scatterName+
" Scatter Residual (Jacobian)");
319 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
324 fieldIds_.resize(scatterFields_.size());
326 const Workset & workset_0 = (*d.worksets_)[0];
327 std::string blockId = this->wda(workset_0).
block_id;
330 for(std::size_t fd=0;fd<scatterFields_.size();++fd) {
332 std::string fieldName = fieldMap_->find(scatterFields_[fd].fieldTag().name())->second;
333 fieldIds_[fd] = globalIndexer_->getFieldNum(fieldName);
335 int fieldNum = fieldIds_[fd];
336 const std::vector<int> &
offsets = globalIndexer_->getGIDFieldOffsets(blockId,fieldNum);
337 scratch_offsets_[fd] = PHX::View<int*>(
"offsets",
offsets.size());
338 Kokkos::deep_copy(scratch_offsets_[fd], Kokkos::View<const int*, Kokkos::HostSpace, Kokkos::MemoryUnmanaged>(
offsets.data(),
offsets.size()));
341 my_derivative_size_ = globalIndexer_->getElementBlockGIDCount(blockId);
342 if (Teuchos::nonnull(workset_0.
other)) {
343 auto otherBlockId = workset_0.
other->block_id;
344 other_derivative_size_ = globalIndexer_->getElementBlockGIDCount(otherBlockId);
346 scratch_lids_ = PHX::View<LO**>(
"lids",scatterFields_[0].extent(0),
347 my_derivative_size_ + other_derivative_size_);
351 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
358 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(d.gedc->getDataObject(globalDataKey_));
360 if(tpetraContainer_==Teuchos::null) {
362 Teuchos::RCP<LinearObjContainer> loc = Teuchos::rcp_dynamic_cast<
LOCPair_GlobalEvaluationData>(d.gedc->getDataObject(globalDataKey_),
true)->getGhostedLOC();
363 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(loc);
372 template <
typename ScalarT,
typename LO,
typename GO,
typename NodeT,
typename LocalMatrixT>
373 class ScatterResidual_Jacobian_Functor {
375 typedef typename PHX::Device execution_space;
376 typedef PHX::MDField<const ScalarT,Cell,NODE>
FieldType;
379 Kokkos::View<double**, Kokkos::LayoutLeft,PHX::Device>
r_data;
387 KOKKOS_INLINE_FUNCTION
388 void operator()(
const unsigned int cell)
const 391 typename Sacado::ScalarType<ScalarT>::type vals[256];
392 int numIds =
lids.extent(1);
394 for(
int i=0;i<numIds;i++)
395 cLIDs[i] =
lids(cell,i);
398 for(std::size_t basis=0; basis <
offsets.extent(0); basis++) {
399 typename FieldType::array_type::reference_type scatterField =
field(cell,basis);
401 LO lid =
lids(cell,offset);
405 Kokkos::atomic_add(&
r_data(lid,0), scatterField.val());
408 for(
int sensIndex=0;sensIndex<numIds;++sensIndex)
409 vals[sensIndex] = scatterField.fastAccessDx(sensIndex);
412 jac.sumIntoValues(lid, cLIDs,numIds, vals,
true,
true);
417 template <
typename ScalarT,
typename LO,
typename GO,
typename NodeT>
418 class ScatterResidual_Residual_Functor {
420 typedef typename PHX::Device execution_space;
421 typedef PHX::MDField<const ScalarT,Cell,NODE>
FieldType;
423 Kokkos::View<double**, Kokkos::LayoutLeft,PHX::Device>
r_data;
425 PHX::View<const LO**>
lids;
430 KOKKOS_INLINE_FUNCTION
431 void operator()(
const unsigned int cell)
const 435 for(std::size_t basis=0; basis <
offsets.extent(0); basis++) {
437 LO lid =
lids(cell,offset);
438 Kokkos::atomic_add(&
r_data(lid,0),
field(cell,basis));
448 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
455 std::string blockId = this->wda(workset).block_id;
457 Teuchos::RCP<typename LOC::VectorType> r = tpetraContainer_->
get_f();
459 globalIndexer_->getElementLIDs(this->wda(workset).cell_local_ids_k,scratch_lids_);
461 ScatterResidual_Residual_Functor<ScalarT,LO,GO,NodeT> functor;
462 functor.r_data = r->getLocalViewDevice(Tpetra::Access::ReadWrite);
463 functor.lids = scratch_lids_;
466 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
467 functor.offsets = scratch_offsets_[fieldIndex];
468 functor.field = scatterFields_[fieldIndex];
470 Kokkos::parallel_for(workset.num_cells,functor);
475 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
481 typedef typename LOC::CrsMatrixType::local_matrix_device_type LocalMatrixT;
484 std::string blockId = this->wda(workset).block_id;
486 Teuchos::RCP<typename LOC::VectorType> r = tpetraContainer_->
get_f();
487 Teuchos::RCP<typename LOC::CrsMatrixType> Jac = tpetraContainer_->get_A();
492 if (Teuchos::nonnull(workset.other)) {
493 auto my_scratch_lids = Kokkos::subview(scratch_lids_,Kokkos::ALL,std::make_pair(0,my_derivative_size_));
494 globalIndexer_->getElementLIDs(this->wda(workset).cell_local_ids_k,my_scratch_lids);
495 auto other_scratch_lids = Kokkos::subview(scratch_lids_,Kokkos::ALL,std::make_pair(my_derivative_size_,my_derivative_size_ + other_derivative_size_));
496 globalIndexer_->getElementLIDs(workset.other->cell_local_ids_k,other_scratch_lids);
499 globalIndexer_->getElementLIDs(this->wda(workset).cell_local_ids_k,scratch_lids_);
502 ScatterResidual_Jacobian_Functor<ScalarT,LO,GO,NodeT,LocalMatrixT> functor;
503 functor.fillResidual = (r!=Teuchos::null);
504 if(functor.fillResidual)
505 functor.r_data = r->getLocalViewDevice(Tpetra::Access::ReadWrite);
506 functor.jac = Jac->getLocalMatrixDevice();
507 functor.lids = scratch_lids_;
510 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
511 functor.offsets = scratch_offsets_[fieldIndex];
512 functor.field = scatterFields_[fieldIndex];
514 Kokkos::parallel_for(workset.num_cells,functor);
panzer::Traits::Tangent::ScalarT ScalarT
FieldType
The type of discretization to use for a field pattern.
Teuchos::RCP< WorksetDetails > other
Kokkos::View< double **, Kokkos::LayoutLeft, PHX::Device > r_data
std::string block_id
DEPRECATED - use: getElementBlock()
Pushes residual values into the residual vector for a Newton-based solve.
PHX::View< const int * > offsets
const Teuchos::RCP< VectorType > get_f() const
PHX::View< const LO ** > lids