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"
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>
74 : globalIndexer_(indexer)
75 , globalDataKey_(
"Residual Scatter Container")
77 std::string scatterName = p.
get<std::string>(
"Scatter Name");
82 const std::vector<std::string>& names =
92 scatterFields_.resize(names.size());
93 scratch_offsets_.resize(names.size());
94 for (std::size_t eq = 0; eq < names.size(); ++eq) {
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());
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) {
149 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(loc);
158 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
162 : globalIndexer_(indexer)
163 , globalDataKey_(
"Residual Scatter Container")
165 std::string scatterName = p.
get<std::string>(
"Scatter Name");
170 const std::vector<std::string>& names =
180 scatterFields_.resize(names.size());
181 for (std::size_t eq = 0; eq < names.size(); ++eq) {
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++) {
229 rcp_dynamic_cast<LOC>(d.gedc->getDataObject(activeParameters[i]),
true)->get_f();
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>
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");
290 const std::vector<std::string>& names =
300 scatterFields_.resize(names.size());
301 scratch_offsets_.resize(names.size());
302 for (std::size_t eq = 0; eq < names.size(); ++eq) {
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());
341 my_derivative_size_ = globalIndexer_->getElementBlockGIDCount(blockId);
343 auto otherBlockId = workset_0.
other->block_id;
344 other_derivative_size_ = globalIndexer_->getElementBlockGIDCount(otherBlockId);
347 "lids", scatterFields_[0].extent(0), my_derivative_size_ + other_derivative_size_ );
349 "vals", scatterFields_[0].extent(0), my_derivative_size_ + other_derivative_size_ );
353 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
360 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(d.gedc->getDataObject(globalDataKey_));
362 if(tpetraContainer_==Teuchos::null) {
365 tpetraContainer_ = Teuchos::rcp_dynamic_cast<LOC>(loc);
374 template <
typename ScalarT,
typename LO,
typename GO,
typename NodeT,
typename LocalMatrixT>
375 class ScatterResidual_Jacobian_Functor {
377 typedef typename PHX::Device execution_space;
390 KOKKOS_INLINE_FUNCTION
391 void operator()(
const unsigned int cell)
const
393 int numIds =
lids.extent(1);
396 for(std::size_t basis=0; basis <
offsets.extent(0); basis++) {
397 typename FieldType::array_type::reference_type scatterField =
field(cell,basis);
399 LO lid =
lids(cell,offset);
403 Kokkos::atomic_add(&
r_data(lid,0), scatterField.val());
406 for(
int sensIndex=0;sensIndex<numIds;++sensIndex)
407 vals(cell,sensIndex) = scatterField.fastAccessDx(sensIndex);
410 jac.sumIntoValues(lid, &
lids(cell,0), numIds, &
vals(cell,0),
true,
true);
415 template <
typename ScalarT,
typename LO,
typename GO,
typename NodeT>
416 class ScatterResidual_Residual_Functor {
418 typedef typename PHX::Device execution_space;
423 PHX::View<const LO**>
lids;
428 KOKKOS_INLINE_FUNCTION
429 void operator()(
const unsigned int cell)
const
433 for(std::size_t basis=0; basis <
offsets.extent(0); basis++) {
435 LO lid =
lids(cell,offset);
436 Kokkos::atomic_add(&
r_data(lid,0),
field(cell,basis));
446 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
453 std::string blockId = this->wda(workset).block_id;
457 globalIndexer_->getElementLIDs(this->wda(workset).cell_local_ids_k,scratch_lids_);
459 ScatterResidual_Residual_Functor<ScalarT,LO,GO,NodeT> functor;
460 functor.r_data = r->getLocalViewDevice(Tpetra::Access::ReadWrite);
461 functor.lids = scratch_lids_;
464 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
465 functor.offsets = scratch_offsets_[fieldIndex];
466 functor.field = scatterFields_[fieldIndex];
468 Kokkos::parallel_for(workset.num_cells,functor);
473 template<
typename TRAITS,
typename LO,
typename GO,
typename NodeT>
479 typedef typename LOC::CrsMatrixType::local_matrix_device_type LocalMatrixT;
482 std::string blockId = this->wda(workset).block_id;
491 auto my_scratch_lids = Kokkos::subview(scratch_lids_,Kokkos::ALL,std::make_pair(0,my_derivative_size_));
492 globalIndexer_->getElementLIDs(this->wda(workset).cell_local_ids_k,my_scratch_lids);
493 auto other_scratch_lids = Kokkos::subview(scratch_lids_,Kokkos::ALL,std::make_pair(my_derivative_size_,my_derivative_size_ + other_derivative_size_));
494 globalIndexer_->getElementLIDs(workset.other->cell_local_ids_k,other_scratch_lids);
497 globalIndexer_->getElementLIDs(this->wda(workset).cell_local_ids_k,scratch_lids_);
500 ScatterResidual_Jacobian_Functor<ScalarT,LO,GO,NodeT,LocalMatrixT> functor;
501 functor.fillResidual = (r!=Teuchos::null);
502 if(functor.fillResidual)
503 functor.r_data = r->getLocalViewDevice(Tpetra::Access::ReadWrite);
504 functor.jac = Jac->getLocalMatrixDevice();
505 functor.lids = scratch_lids_;
506 functor.vals = scratch_vals_;
509 for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
510 functor.offsets = scratch_offsets_[fieldIndex];
511 functor.field = scatterFields_[fieldIndex];
513 Kokkos::parallel_for(workset.num_cells,functor);
panzer::Traits::Tangent::ScalarT ScalarT
T & get(const std::string &name, T def_value)
bool nonnull(const std::shared_ptr< T > &p)
FieldType
The type of discretization to use for a field pattern.
PHX::View< const int * > offsets
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos::RCP< WorksetDetails > other
Kokkos::View< typename Sacado::ScalarType< ScalarT >::type **, Kokkos::LayoutRight, PHX::Device > vals
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::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
Kokkos::View< const LO **, Kokkos::LayoutRight, PHX::Device > lids